Packages & Versions Used

# Table of packages
kable(table[-1,], format = "html", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Package Title Maintainer Version URL
foreign Read Data Stored by ‘Minitab’, ‘S’, ‘SAS’, ‘SPSS’, ‘Stata’, ‘Systat’, ‘Weka’, ‘dBase’, … R Core Team <; 0.8-80 NA
tidyverse Easily Install and Load the ‘Tidyverse’ Hadley Wickham <; 1.3.0 http://tidyverse.tidyverse.org,
MplusAutomation An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus Michael Hallquist <; 0.7-3 NA
interactions Comprehensive, User-Friendly Toolkit for Probing Interactions Jacob A. Long <; 1.1.3 NA
ggpubr ‘ggplot2’ Based Publication Ready Plots Alboukadel Kassambara <; 0.4.0 NA
psych Procedures for Psychological, Psychometric, and Personality Research William Revelle <; 2.0.7 NA
table1 Tables of Descriptive Statistics in HTML Benjamin Rich <; 1.2 NA
MatchIt Nonparametric Preprocessing for Parametric Causal Inference Kosuke Imai <; 3.0.2 NA
apaTables Create American Psychological Association (APA) Style Tables David Stanley <; 2.0.5 NA
car Companion to Applied Regression John Fox <; 3.0-9 https://r-forge.r-project.org/projects/car/,
MASS Support Functions and Datasets for Venables and Ripley’s MASS Brian Ripley <; 7.3-51.6 NA
data.table Extension of data.frame Matt Dowle <; 1.13.0 http://r-datatable.com, https://Rdatatable.gitlab.io/data.table,
kableExtra Construct Complex Table with ‘kable’ and Pipe Syntax Hao Zhu <; 1.1.0 http://haozhu233.github.io/kableExtra/,

Import Data

FullData = read.spss('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/FFCW_core_ffid_transfer.sav',to.data.frame=TRUE, use.value.labels = FALSE)

Preparing the Data

The composite scores for violence exposure and social deprivation are calculated exactly how they’re done in Hein et al., 2020 in SCAN.

Creating Some Initial Flags

###Creating a variable to indicate whether child lived with mother most of the time at each time point

FullData$age3cust <- ifelse(FullData$m3a2 == 1 |FullData$m3a2 == 2, "1", "0")
FullData$age5cust <- ifelse(FullData$m4a2 == 1 |FullData$m4a2 == 2, "1", "0")
FullData$age9cust <- ifelse(FullData$m5a2 == 1 |FullData$m5a2 == 2, "1", "0")

###Creating a variable to indicate whether mom was currently married, living with, or romantically involved with baby’s father at age 9 (since m5c4 doesn’t seem to be recognized) ####FFCWS documentation indicates that they would assign a 1 / yes for the following responses to a4 - 1(married), 4(cohabiting or living together), 5(romantically involved but living apart), -1(refused), -2(don’t know) ####Based on conversations with Colter Mitchell about refusing itself being of interest, will instead only include1, 4, or 5

FullData$con.m5c4 = ifelse(FullData$m5a4 == 1 | FullData$m5a4 == 4 | FullData$m5a4 == 5, "1", "0")

###Creating a variable to indicate whether mom was currently in any relationship at each time point

FullData$age3rel = ifelse(FullData$m3d6 == 1 |FullData$m3e2 == 1, "1", "0")
FullData$age5rel = ifelse(FullData$m4d5 == 1 |FullData$m4e2 == 1, "1", "0")
FullData$age9rel = ifelse(FullData$con.m5c4 == 1 |FullData$m5d2 == 1, "1", "0")

###Creating a flag to indicate that the mom said she was not in a relationship - with BF or CP

FullData$age3.norel = ifelse(FullData$m3d6 == 2 & FullData$m3e2 == 2, "1", "0")
FullData$age5.norel = ifelse(FullData$m4d5 == 2 & FullData$m4e2 == 2, "1", "0")
FullData$age9.norel = ifelse(FullData$con.m5c4 == 2 & FullData$m5d2 == 2, "1", "0")

###Creating a flag to indicate that the mom refused to answer question about relationship with BF

FullData$age3.refusebf = ifelse(FullData$m3a4 == -1, "1", "0")
FullData$age5.refusebf = ifelse(FullData$m4a4 == -1, "1", "0")
FullData$age9.refusebf = ifelse(FullData$m5a4 == -1, "1", "0")

###Creating a flag to indicate that the mom refused to answer question about relationship with a CP

FullData$age3.refusecp = ifelse(FullData$m3e2 == -1, "1", "0")
FullData$age5.refusecp = ifelse(FullData$m4e2 == -1, "1", "0")
FullData$age9.refusecp = ifelse(FullData$m5d2 == -1, "1", "0")

Recode for Missing

-9 (not in wave), -6 (skip), -3(missing), -2(don’t know), -1 (refuse)

is.na(FullData) <- FullData == -9
is.na(FullData) <- FullData == -8
is.na(FullData) <- FullData == -6
is.na(FullData) <- FullData == -5
is.na(FullData) <- FullData == -3
is.na(FullData) <- FullData == -2
is.na(FullData) <- FullData == -1

Select Variables

Violence Exposure and Social Deprivation Composite Scores The items used in the violence exposure and social deprivation composites come from Hein et al., 2020, SCAN are are the same as those in Goetschius et al., 2020, DCN, Goetschius et al., 2020, JAMA Network Open, and Peckins et al., 2019, PNE. I haven’t repeated where they come from here, but their origins can be found in those papers. Additionally, specific papers that T.Hein selected items based on should be referenced in the code below.

School Connectedness 9 & 15 Items for this come from the Panel Study on Income Dynamics Child Development Supplement (PSID-CDS) (The Panel Study of Income Dynamics Child Development Supplement: User Guide for CDS-III, 2010)

Child Internalizing and Externalizing Symptoms Items for both of these came from the Child Behavior Checklist (Achenbach & Edelbrock, 1983). The internalizing items came from the anxious/depressed, withdrawn/depressed, and somatic symptoms subscales and the externalizing items came from the rule-breaking and aggressive behavior subscales.

Positive Adolescent Function Items from this come from the EPOCH Measure of Adolescent Wellbeing (Kern et al., 2016)

Adolescent Internalizing Symptoms Items for this come from: - Center for Epidemiologic Studies Depression Scale (CES-D) (Radloff, 1977) - Brief Symptom Inventory 18 (BSI-18) (Derogatis & Savitz, 2000)

Adolescent Externalizing Symptoms Items for this come from: - Delinquency: National Longitudinal Study of Adolescent Health (Add Health - (Harris, 2013) - Dickman’s Impulsivity scale (Dickman, 1990) - Binary substance use variables (alcohol-2+ drinks, cigarettes, marijuana, illegal drugs other than marijuana, prescription drugs w/out prescription)

Covariates Most of the covariates are pretty self-explanatory, but one requires a bit more information. To get at income I used the poverty ratio that’s constructed in the Fragile Families dataset. Because we use data from most of the waves and I don’t want to have multiple covariates representing income in my SEM for model degrees of freedom purposes, I calculate the average poverty ratio and use that as the covariate. According to the Fragile Families website, the poverty ratio is the ratio of total household income, as defined in c_1hhinc, to the official poverty thresholds, designated by the U.S. Census Bureau. The thresholds in c_1povca vary by family composition and year. At each wave, we used the poverty thresholds for the year preceding the interview. We calculated separate thresholds based on mother and father reports of household size and composition. However, calculations for married/cohabiting mothers and fathers rely on mother reports of household size and composition. A small number of missing values (don’t know, refused) were treated as 0 in household membership counts.Also included is m1city.

Primary Caregiver Relationship Used this to determine who the primary caregiver generally was for the reported data at the study waves included in the violence exposure and social deprivation composites as well as the externalizing and internalizing at age 9.

Data.ThreatDepVars = FullData %>% 
  dplyr::select( m1b2, cm1edu,cm2md_case_con,cm2md_case_lib,m2b17a,m2b17b,m2b17c,m2b17d,m2b17e,m2b17f,IH3j10,IH3j11,IH3j13,IH3j14,IH3j15,IH3j16,IH3j17,IH3j18,IH3j19,IH3j3,IH3j4,IH3j6,IH3j7,IH3j8,IH3j9,IH3k1a,IH3k1b,IH3k1c,IH3k1d,IH3k1e,IH3k2a,IH3k2b,IH3k2c,IH3k2d,IH3k2e,IH3l1,IH3l2,IH3l3,IH3l4,IH3l5,IH3l6,IH3l7,m3a2,m3a4,m3d6,m3d9n1,m3d9p,m3d9p1,m3d10,m3d10d,m3d1f,m3d7a,m3d7b,m3d7c,m3d7d,m3d7e,m3d7f,m3d7g,m3d7h,m3d7i,m3d7j,m3d7k,m3d7l,m3d7n,m3d7n1,m3d7p,m3d7p1,m3e2,m3e21b,m3e21f,m3e23a,m3e23b,m3e23c,m3e23d,m3e23e,m3e23f,m3e23g,m3e23h,m3e23i,m3e23j,m3e23k,m3e23l,m3e24,m3h3,m3h3a,m3h4,m3h5,m3h6,m3i6c, m3i6f,m3i23f,m3j36j,m3j36k,m3k26a, m3k26b,IH3r10,IH3s2,IH3s3,IH3s4,IH5g3,IH5g4,IH5g6,IH5g7,IH5g8,IH5g9,IH5g10,IH5g11,IH5g13,IH5g14,IH5g15,IH5g16,IH5g17,IH5g18,IH5g19,IH5h1,IH5h2,IH5h3,IH5h4,IH5h5,IH5h6,IH5h7,m4a2,m4a4,m4d1b,m4d1f,m4d5,m4d7a,m4d7b,m4d7c,m4d7d,m4d7e,m4d7f,m4d7g,m4d7h,m4d7i,m4d7j,m4d7o,m4d7p,m4d10,m4d10a,m4d10e,m4e2,m4e21b,m4e21f,m4e23a,m4e23b,m4e23c,m4e23d,m4e23e,m4e23f,m4e23g,m4e23h,m4e23i,m4e23j,m4e23o,m4e23p,m4e23q,m4e24,m4e24d,m4h3,m4h3a,m4h4,m4h5,m4h6,m4i0m1,m4i0m2,m4i0m3,m4i0m4,m4i0m5,m4i0n1,m4i0n2,m4i0n3,m4i0n4,m4i0n5,m4i23g,m4i23i,m4j22j,m4j22k,m4k25a,m4k25b,IH5r10,IH5s2,IH5s3,IH5s4,IH5s5,o5c10,o5d2,o5d3,o5d4,o5d5,n5g1f,n5g1h,hv5_dsae,hv5_dspr,hv5_dsraw,hv5_dsss,hv5_ppvtae,hv5_ppvtpr,hv5_ppvtraw,hv5_ppvtss,p5m2a,p5m2b,p5m2c,p5m2d,p5m2e,p5m3a,p5m3b,p5m3c,p5m3d,p5m3e,p5m5a,p5m5b,p5m5c,m5a2,m5a4,m5c2b,m5c2f,m5c6a,m5c6b,m5c6c,m5c6d,m5c6e,m5c6f,m5c6g,m5c6h,m5c6i,m5c6j,m5c6o,m5c6p,m5c7,m5c8a,m5c8e,m5d2,m5d19b,m5d19f,m5d20a,m5d20b,m5d20c,m5d20d,m5d20e,m5d20f,m5d20g,m5d20h,m5d20i,m5d20j,m5d20o,m5d20p,m5d21,m5d22,m5d22d,m5e3,m5e3a,m5e4,m5e5,m5e6,m5g21a,m5g21b,m5g21c,m5g21d,m5g21e,m5g21f,m5g21g,m5g21h,m5g21i,m5g21k,m5i25a,m5i25b,p5q1c,p5q1d,p5q1f,p5q1g,p5q1h,p5q1i,p5q1j,p5q1k,p5q1m,p5q1n,p5q2a,p5q2b,p5q2c,p5q2d,p5q2e,ff_id,m1city,age3cust:age9.refusecp)

Data.PAF = FullData %>% 
  dplyr::select("ff_id", "k6d2b", "k6d2e", "k6d2f", "k6d2g", "k6d2h", "k6d2i", "k6d2k", "k6d2l", "k6d2m", "k6d2o", "k6d2s", "k6d2u", "k6d2v", "k6d2w", "k6d2y", "k6d2aa", "k6d2ad", "k6d2ae", "k6d2af", "k6d2ah")

Data.AnxDep = FullData %>% 
  dplyr::select(ff_id, k6d2ag, k6d2ai, k6d2d, k6d2j, k6d2t, k6d2ac, k6d2ak, k6d2c, k6d2n, k6d2x, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66)

Data.SchoolConn = FullData %>% 
  dplyr::select("ff_id", "k5e1a", "k5e1b", "k5e1c", "k5e1d", "k6b1a", "k6b1b", "k6b1c", "k6b1d")

Extern_Variables = FullData %>% 
  dplyr::select(ff_id, k6d61a:k6d61m, k6d2a, k6d2p, k6d2r, k6d2z, k6d2ab, k6d2aj, k6d40, k6d48, k6f63, k6f68, k6f74, p6b35, p6b37:p6b39, p6b41:p6b45, p6b57, p6b59, p6b49:p6b51, p6b60:p6b64, p6b67)

CBCL = FullData %>% 
  dplyr::select(ff_id,p5q3m,p5q3ab,p5q3ac,p5q3ad,p5q3ae,p5q3af,p5q3ah,p5q3ar,p5q3av,p5q3ax,p5q3bq,p5q3ck,p5q3db,p5q3e,p5q3ao,p5q3bk,p5q3bo,p5q3bu,p5q3cu,p5q3cv,p5q3da,p5q3as,p5q3au,p5q3aw,p5q3az,p5q3bb1,p5q3bb2,p5q3bb3,p5q3bb4,p5q3bb5,p5q3bb6,p5q3bb7,p5q3b,p5q3x,p5q3aa,p5q3al,p5q3ap,p5q3bi,p5q3bm,p5q3br,p5q3bs,p5q3bz,p5q3ca,p5q3cj,p5q3cp,p5q3cr,p5q3ct,p5q3cx,p5q3cy,p5q3c,p5q3o,p5q3r,p5q3s,p5q3t,p5q3u,p5q3v,p5q3aj,p5q3bc,p5q3bn,p5q3cf,p5q3cg,p5q3ch,p5q3ci,p5q3cn,p5q3co,p5q3cq,p5q3cw)

Demo_Variables = FullData %>% 
  dplyr::select(ff_id, ck6ethrace, cm1bsex, cm1inpov, cm2povco, cm3povco, cm4povco, cm5povco, cp6povco, m1city)

PC_Info = FullData %>% 
  dplyr::select(ff_id, IH3int5, IH5int5, pcg5idstat)

CBCL_15 = FullData %>% 
  dplyr::select(ff_id, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66, p6b35, p6b37, p6b38, p6b39, p6b41, p6b42, p6b43, p6b44, p6b45, p6b57, p6b59, p6b49, p6b50, p6b51, p6b60, p6b61, p6b62, p6b63, p6b64, p6b67)

Recode Variables

School Connectedness

We need an extra NA code for the school connectedness variables (particularly the age 15 variables) – 7 on the school connectedness scale means (NA - Homeschooled)

is.na(Data.SchoolConn) <- Data.SchoolConn == 7

Data.SchoolConn$k6b1a_r = ((Data.SchoolConn$k6b1a - 5)*-1)
Data.SchoolConn$k6b1b_r = ((Data.SchoolConn$k6b1b - 5)*-1)
Data.SchoolConn$k6b1c_r = ((Data.SchoolConn$k6b1c - 5)*-1)
Data.SchoolConn$k6b1d_r = ((Data.SchoolConn$k6b1d - 5)*-1)

Teen Report Internalizing

Data.AnxDep$k6d2ag_r = ((Data.AnxDep$k6d2ag-5)*-1)
Data.AnxDep$k6d2ai_r = ((Data.AnxDep$k6d2ai-5)*-1)
Data.AnxDep$k6d2d_r = ((Data.AnxDep$k6d2d-5)*-1)
Data.AnxDep$k6d2j_r = ((Data.AnxDep$k6d2j-5)*-1)
Data.AnxDep$k6d2t_r = ((Data.AnxDep$k6d2t-5)*-1)
Data.AnxDep$k6d2ac_r = ((Data.AnxDep$k6d2ac-5)*-1)
Data.AnxDep$k6d2ak_r = ((Data.AnxDep$k6d2ak-5)*-1)
Data.AnxDep$k6d2c_r = ((Data.AnxDep$k6d2c-5)*-1)
Data.AnxDep$k6d2n_r = ((Data.AnxDep$k6d2n-5)*-1)
Data.AnxDep$k6d2x_r = ((Data.AnxDep$k6d2x-5)*-1)

Teen Report Externalizing

Extern_Variables = Extern_Variables %>% 
  mutate(k6d2a_r = (k6d2a-5)*-1,
         k6d2p_r = (k6d2p-5)*-1,
         k6d2r_r = (k6d2r-5)*-1,
         k6d2z_r = (k6d2z-5)*-1,
         k6d2ab_r = (k6d2ab-5)*-1,
         k6d2aj_r = (k6d2aj-5)*-1,
         k6d40_r = (k6d40-3)*-1,
         k6d48_r = (k6d48-3)*-1,
         k6f63_r = (k6f63-3)*-1,
         k6f68_r = (k6f68-3)*-1,
         k6f74_r = (k6f74-2)*-1)

Positive Adolescent Function

Data.PAF$k6d2b_r = ((Data.PAF$k6d2b-5)*-1)
Data.PAF$k6d2b_r = 5-Data.PAF$k6d2b
Data.PAF$k6d2e_r = ((Data.PAF$k6d2e-5)*-1)
Data.PAF$k6d2f_r = ((Data.PAF$k6d2f-5)*-1)
Data.PAF$k6d2g_r = ((Data.PAF$k6d2g-5)*-1)
Data.PAF$k6d2h_r = ((Data.PAF$k6d2h-5)*-1)
Data.PAF$k6d2i_r = ((Data.PAF$k6d2i-5)*-1)
Data.PAF$k6d2k_r = ((Data.PAF$k6d2k-5)*-1)
Data.PAF$k6d2l_r = ((Data.PAF$k6d2l-5)*-1)
Data.PAF$k6d2m_r = ((Data.PAF$k6d2m-5)*-1)
Data.PAF$k6d2o_r = ((Data.PAF$k6d2o-5)*-1)
Data.PAF$k6d2s_r = ((Data.PAF$k6d2s-5)*-1)
Data.PAF$k6d2u_r = ((Data.PAF$k6d2u-5)*-1)
Data.PAF$k6d2v_r = ((Data.PAF$k6d2v-5)*-1)
Data.PAF$k6d2w_r = ((Data.PAF$k6d2w-5)*-1)
Data.PAF$k6d2y_r = ((Data.PAF$k6d2y-5)*-1)
Data.PAF$k6d2aa_r = ((Data.PAF$k6d2aa-5)*-1)
Data.PAF$k6d2ad_r = ((Data.PAF$k6d2ad-5)*-1)
Data.PAF$k6d2ae_r = ((Data.PAF$k6d2ae-5)*-1)
Data.PAF$k6d2af_r = ((Data.PAF$k6d2af-5)*-1)
Data.PAF$k6d2ah_r = ((Data.PAF$k6d2ah-5)*-1)

Calculate Violence Exposure and Social Dep Composites

##Abuse ###Item selection based on Hunt, Slack, and Berger 2016 ###Convert all items with 7 (yes but not in last year) and 8 (this has never happened) to 0

Data.ThreatDepVars$IH3j7<-car::recode(Data.ThreatDepVars$IH3j7 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j4<-car::recode(Data.ThreatDepVars$IH3j4 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j11<-car::recode(Data.ThreatDepVars$IH3j11 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j13<-car::recode(Data.ThreatDepVars$IH3j13 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j3<-car::recode(Data.ThreatDepVars$IH3j3 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j14<-car::recode(Data.ThreatDepVars$IH3j14 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j9<-car::recode(Data.ThreatDepVars$IH3j9 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j8<-car::recode(Data.ThreatDepVars$IH3j8 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j6<-car::recode(Data.ThreatDepVars$IH3j6 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH3j10<-car::recode(Data.ThreatDepVars$IH3j10 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g7<-car::recode(Data.ThreatDepVars$IH5g7 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g4<-car::recode(Data.ThreatDepVars$IH5g4 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g11<-car::recode(Data.ThreatDepVars$IH5g11 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g13<-car::recode(Data.ThreatDepVars$IH5g13 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g3<-car::recode(Data.ThreatDepVars$IH5g3 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g14<-car::recode(Data.ThreatDepVars$IH5g14 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g9<-car::recode(Data.ThreatDepVars$IH5g9 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g8<-car::recode(Data.ThreatDepVars$IH5g8 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g6<-car::recode(Data.ThreatDepVars$IH5g6 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$IH5g10<-car::recode(Data.ThreatDepVars$IH5g10 , "c(7,8)=0", as.numeric=TRUE,as.factor=FALSE)
Data.ThreatDepVars$p5q1g<-car::recode(Data.ThreatDepVars$p5q1g, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1d<-car::recode(Data.ThreatDepVars$p5q1d, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1k<-car::recode(Data.ThreatDepVars$p5q1k, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1m<-car::recode(Data.ThreatDepVars$p5q1m, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1c<-car::recode(Data.ThreatDepVars$p5q1c, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1n<-car::recode(Data.ThreatDepVars$p5q1n, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1i<-car::recode(Data.ThreatDepVars$p5q1i, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1h<-car::recode(Data.ThreatDepVars$p5q1h, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1f<-car::recode(Data.ThreatDepVars$p5q1f, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)
Data.ThreatDepVars$p5q1j<-car::recode(Data.ThreatDepVars$p5q1j, "c(7,8) = 0",  as.numeric = TRUE, as.factor =FALSE)

###Calculating physical abuse, emotional abuse, and abuse scores

Data.ThreatDepVars$Age3.PhysicalAbuse = Data.ThreatDepVars$IH3j7 + Data.ThreatDepVars$IH3j4 + Data.ThreatDepVars$IH3j11 + Data.ThreatDepVars$IH3j13 + Data.ThreatDepVars$IH3j3
Data.ThreatDepVars$Age3.EmotionalAbuse =Data.ThreatDepVars$IH3j14 + Data.ThreatDepVars$IH3j9 + Data.ThreatDepVars$IH3j8 + Data.ThreatDepVars$IH3j6 + Data.ThreatDepVars$IH3j10
Data.ThreatDepVars$Age3.Abuse = Data.ThreatDepVars$Age3.PhysicalAbuse + Data.ThreatDepVars$Age3.EmotionalAbuse
Data.ThreatDepVars$Age5.PhysicalAbuse = Data.ThreatDepVars$IH5g7 + Data.ThreatDepVars$IH5g4 + Data.ThreatDepVars$IH5g11 + Data.ThreatDepVars$IH5g13 + Data.ThreatDepVars$IH5g3
Data.ThreatDepVars$Age5.EmotionalAbuse = Data.ThreatDepVars$IH5g14 + Data.ThreatDepVars$IH5g9 + Data.ThreatDepVars$IH5g8 + Data.ThreatDepVars$IH5g6 + Data.ThreatDepVars$IH5g10
Data.ThreatDepVars$Age5.Abuse = Data.ThreatDepVars$Age5.PhysicalAbuse + Data.ThreatDepVars$Age5.EmotionalAbuse
Data.ThreatDepVars$Age9.PhysicalAbuse = Data.ThreatDepVars$p5q1g + Data.ThreatDepVars$p5q1d + Data.ThreatDepVars$p5q1k + Data.ThreatDepVars$p5q1m + Data.ThreatDepVars$p5q1c
Data.ThreatDepVars$Age9.EmotionalAbuse = Data.ThreatDepVars$p5q1n + Data.ThreatDepVars$p5q1i + Data.ThreatDepVars$p5q1h + Data.ThreatDepVars$p5q1f + Data.ThreatDepVars$p5q1j
Data.ThreatDepVars$Age9.Abuse = Data.ThreatDepVars$Age9.PhysicalAbuse + Data.ThreatDepVars$Age9.EmotionalAbuse

Internal Reliability

I have this in here because I’ve been asked previously about how the items relate to each other and about their alpha scores.

# Age 3 Abuse
Data.ThreatDepVars %>% 
  dplyr::select("IH3j7", "IH3j4", "IH3j11", "IH3j13", "IH3j3", "IH3j14", "IH3j9", "IH3j8", "IH3j6", "IH3j10") %>% 
  psych::alpha(.)

# Age 5 Abuse

Data.ThreatDepVars %>%
  dplyr::select("IH5g7", "IH5g4", "IH5g11", "IH5g13", "IH5g3", "IH5g14", "IH5g9", "IH5g8", "IH5g6", "IH5g10") %>% 
  psych::alpha(.)

# Age 9 Abuse

Data.ThreatDepVars %>% 
  dplyr::select("p5q1g", "p5q1d", "p5q1k", "p5q1m", "p5q1c", "p5q1n", "p5q1i", "p5q1h", "p5q1f", "p5q1j") %>% 
  psych::alpha(.)

##IPV ###Item selection based on Hunt, Slack, & Berger 2016 ###Recode items - 3 (never) becomes 0, 2 becomes 1, 1 becomes 2 - this matches abuse, neglect, and CV where 0 = never

Data.ThreatDepVars$m3d7h<-car::recode(Data.ThreatDepVars$m3d7h, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3d7i<-car::recode(Data.ThreatDepVars$m3d7i, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3d7e<-car::recode(Data.ThreatDepVars$m3d7e, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3d7f<-car::recode(Data.ThreatDepVars$m3d7f, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3d7g<-car::recode(Data.ThreatDepVars$m3d7g, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3d7j<-car::recode(Data.ThreatDepVars$m3d7j, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7h<-car::recode(Data.ThreatDepVars$m4d7h, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7i<-car::recode(Data.ThreatDepVars$m4d7i, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7e<-car::recode(Data.ThreatDepVars$m4d7e, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7f<-car::recode(Data.ThreatDepVars$m4d7f, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7g<-car::recode(Data.ThreatDepVars$m4d7g, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4d7j<-car::recode(Data.ThreatDepVars$m4d7j, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6h<-car::recode(Data.ThreatDepVars$m5c6h, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6i<-car::recode(Data.ThreatDepVars$m5c6i, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6e<-car::recode(Data.ThreatDepVars$m5c6e, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6f<-car::recode(Data.ThreatDepVars$m5c6f, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6g<-car::recode(Data.ThreatDepVars$m5c6g, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5c6j<-car::recode(Data.ThreatDepVars$m5c6j, "3 = 0 ; 2 = 1 ; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23h<-car::recode(Data.ThreatDepVars$m3e23h, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23i<-car::recode(Data.ThreatDepVars$m3e23i, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23e<-car::recode(Data.ThreatDepVars$m3e23e, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23f<-car::recode(Data.ThreatDepVars$m3e23f, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23g<-car::recode(Data.ThreatDepVars$m3e23g, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m3e23j<-car::recode(Data.ThreatDepVars$m3e23j, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23h<-car::recode(Data.ThreatDepVars$m4e23h, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23i<-car::recode(Data.ThreatDepVars$m4e23i, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23e<-car::recode(Data.ThreatDepVars$m4e23e, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23f<-car::recode(Data.ThreatDepVars$m4e23f, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23g<-car::recode(Data.ThreatDepVars$m4e23g, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m4e23j<-car::recode(Data.ThreatDepVars$m4e23j, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20h<-car::recode(Data.ThreatDepVars$m5d20h, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20i<-car::recode(Data.ThreatDepVars$m5d20i, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20e<-car::recode(Data.ThreatDepVars$m5d20e, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20f<-car::recode(Data.ThreatDepVars$m5d20f, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20g<-car::recode(Data.ThreatDepVars$m5d20g, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE) 
Data.ThreatDepVars$m5d20j<-car::recode(Data.ThreatDepVars$m5d20j, "3 = 0; 2 = 1; 1 = 2", as.numeric = TRUE)

###Calculating physical, emotional, sexual, and total IPV for baby’s father, current partner

Data.ThreatDepVars$BF.IPV.3.Physical = Data.ThreatDepVars$m3d7h + Data.ThreatDepVars$m3d7i
Data.ThreatDepVars$BF.IPV.3.Emotional = Data.ThreatDepVars$m3d7e + Data.ThreatDepVars$m3d7f + Data.ThreatDepVars$m3d7g
Data.ThreatDepVars$BF.IPV.3.Sexual = Data.ThreatDepVars$m3d7j
Data.ThreatDepVars$BF.IPV.3.Total = Data.ThreatDepVars$BF.IPV.3.Physical + Data.ThreatDepVars$BF.IPV.3.Emotional + Data.ThreatDepVars$BF.IPV.3.Sexual
Data.ThreatDepVars$BF.IPV.5.Physical = Data.ThreatDepVars$m4d7h + Data.ThreatDepVars$m4d7i
Data.ThreatDepVars$BF.IPV.5.Emotional = Data.ThreatDepVars$m4d7e + Data.ThreatDepVars$m4d7f + Data.ThreatDepVars$m4d7g
Data.ThreatDepVars$BF.IPV.5.Sexual = Data.ThreatDepVars$m4d7j
Data.ThreatDepVars$BF.IPV.5.Total = Data.ThreatDepVars$BF.IPV.5.Physical + Data.ThreatDepVars$BF.IPV.5.Emotional + Data.ThreatDepVars$BF.IPV.5.Sexual
Data.ThreatDepVars$BF.IPV.9.Physical = Data.ThreatDepVars$m5c6h + Data.ThreatDepVars$m5c6i
Data.ThreatDepVars$BF.IPV.9.Emotional = Data.ThreatDepVars$m5c6e + Data.ThreatDepVars$m5c6f + Data.ThreatDepVars$m5c6g
Data.ThreatDepVars$BF.IPV.9.Sexual = Data.ThreatDepVars$m5c6j
Data.ThreatDepVars$BF.IPV.9.Total = Data.ThreatDepVars$BF.IPV.9.Physical + Data.ThreatDepVars$BF.IPV.9.Emotional + Data.ThreatDepVars$BF.IPV.9.Sexual
Data.ThreatDepVars$CP.IPV.3.Physical = Data.ThreatDepVars$m3e23h + Data.ThreatDepVars$m3e23i
Data.ThreatDepVars$CP.IPV.3.Emotional = Data.ThreatDepVars$m3e23e + Data.ThreatDepVars$m3e23f + Data.ThreatDepVars$m3e23g
Data.ThreatDepVars$CP.IPV.3.Sexual = Data.ThreatDepVars$m3e23j
Data.ThreatDepVars$CP.IPV.3.Total = Data.ThreatDepVars$CP.IPV.3.Physical + Data.ThreatDepVars$CP.IPV.3.Emotional + Data.ThreatDepVars$CP.IPV.3.Sexual
Data.ThreatDepVars$CP.IPV.5.Physical = Data.ThreatDepVars$m4e23h + Data.ThreatDepVars$m4e23i
Data.ThreatDepVars$CP.IPV.5.Emotional = Data.ThreatDepVars$m4e23e + Data.ThreatDepVars$m4e23f + Data.ThreatDepVars$m4e23g
Data.ThreatDepVars$CP.IPV.5.Sexual = Data.ThreatDepVars$m4e23j
Data.ThreatDepVars$CP.IPV.5.Total = Data.ThreatDepVars$CP.IPV.5.Physical + Data.ThreatDepVars$CP.IPV.5.Emotional + Data.ThreatDepVars$CP.IPV.5.Sexual
Data.ThreatDepVars$CP.IPV.9.Physical = Data.ThreatDepVars$m5d20h + Data.ThreatDepVars$m5d20i
Data.ThreatDepVars$CP.IPV.9.Emotional = Data.ThreatDepVars$m5d20e + Data.ThreatDepVars$m5d20f + Data.ThreatDepVars$m5d20g
Data.ThreatDepVars$CP.IPV.9.Sexual = Data.ThreatDepVars$m5d20j
Data.ThreatDepVars$CP.IPV.9.Total = Data.ThreatDepVars$CP.IPV.9.Physical + Data.ThreatDepVars$CP.IPV.9.Emotional + Data.ThreatDepVars$CP.IPV.9.Sexual

###Merge baby’s father, current partner columns

Data.ThreatDepVars$IPV.3.Physical = Data.ThreatDepVars$BF.IPV.3.Physical
Data.ThreatDepVars$IPV.3.Physical[!is.na(Data.ThreatDepVars$CP.IPV.3.Physical)] = Data.ThreatDepVars$CP.IPV.3.Physical[!is.na(Data.ThreatDepVars$CP.IPV.3.Physical)]
Data.ThreatDepVars$IPV.3.Emotional = Data.ThreatDepVars$BF.IPV.3.Emotional
Data.ThreatDepVars$IPV.3.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.3.Emotional)] = Data.ThreatDepVars$CP.IPV.3.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.3.Emotional)]
Data.ThreatDepVars$IPV.3.Sexual = Data.ThreatDepVars$BF.IPV.3.Sexual
Data.ThreatDepVars$IPV.3.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.3.Sexual)] = Data.ThreatDepVars$CP.IPV.3.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.3.Sexual)]
Data.ThreatDepVars$IPV.3.Total = Data.ThreatDepVars$BF.IPV.3.Total
Data.ThreatDepVars$IPV.3.Total[!is.na(Data.ThreatDepVars$CP.IPV.3.Total)] = Data.ThreatDepVars$CP.IPV.3.Total[!is.na(Data.ThreatDepVars$CP.IPV.3.Total)]
Data.ThreatDepVars$IPV.5.Physical = Data.ThreatDepVars$BF.IPV.5.Physical
Data.ThreatDepVars$IPV.5.Physical[!is.na(Data.ThreatDepVars$CP.IPV.5.Physical)] = Data.ThreatDepVars$CP.IPV.5.Physical[!is.na(Data.ThreatDepVars$CP.IPV.5.Physical)]
Data.ThreatDepVars$IPV.5.Emotional = Data.ThreatDepVars$BF.IPV.5.Emotional
Data.ThreatDepVars$IPV.5.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.5.Emotional)] = Data.ThreatDepVars$CP.IPV.5.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.5.Emotional)]
Data.ThreatDepVars$IPV.5.Sexual = Data.ThreatDepVars$BF.IPV.5.Sexual
Data.ThreatDepVars$IPV.5.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.5.Sexual)] = Data.ThreatDepVars$CP.IPV.5.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.5.Sexual)]
Data.ThreatDepVars$IPV.5.Total = Data.ThreatDepVars$BF.IPV.5.Total
Data.ThreatDepVars$IPV.5.Total[!is.na(Data.ThreatDepVars$CP.IPV.5.Total)] = Data.ThreatDepVars$CP.IPV.5.Total[!is.na(Data.ThreatDepVars$CP.IPV.5.Total)]
Data.ThreatDepVars$IPV.9.Physical = Data.ThreatDepVars$BF.IPV.9.Physical
Data.ThreatDepVars$IPV.9.Physical[!is.na(Data.ThreatDepVars$CP.IPV.9.Physical)] = Data.ThreatDepVars$CP.IPV.9.Physical[!is.na(Data.ThreatDepVars$CP.IPV.9.Physical)]
Data.ThreatDepVars$IPV.9.Emotional = Data.ThreatDepVars$BF.IPV.9.Emotional
Data.ThreatDepVars$IPV.9.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.9.Emotional)] = Data.ThreatDepVars$CP.IPV.9.Emotional[!is.na(Data.ThreatDepVars$CP.IPV.9.Emotional)]
Data.ThreatDepVars$IPV.9.Sexual = Data.ThreatDepVars$BF.IPV.9.Sexual
Data.ThreatDepVars$IPV.9.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.9.Sexual)] = Data.ThreatDepVars$CP.IPV.9.Sexual[!is.na(Data.ThreatDepVars$CP.IPV.9.Sexual)]
Data.ThreatDepVars$IPV.9.Total = Data.ThreatDepVars$BF.IPV.9.Total
Data.ThreatDepVars$IPV.9.Total[!is.na(Data.ThreatDepVars$CP.IPV.9.Total)] = Data.ThreatDepVars$CP.IPV.9.Total[!is.na(Data.ThreatDepVars$CP.IPV.9.Total)]

##Community Violence ###Item selection based on Zhang & Anderson 2010 but note they only did age 3 and age 9 does not have victimization items ###Calculate sum scores for community violence witnessing and victimization ###Going to create a sum of what is available at each age to maximize threat index for use with the composite score

Data.ThreatDepVars$CV.Age3.Witness.KillIncluded = Data.ThreatDepVars$IH3l1 + Data.ThreatDepVars$IH3l3 + Data.ThreatDepVars$IH3l5 + Data.ThreatDepVars$IH3l7
Data.ThreatDepVars$CV.Age3.Witness = Data.ThreatDepVars$IH3l1 + Data.ThreatDepVars$IH3l3 + Data.ThreatDepVars$IH3l5
Data.ThreatDepVars$CV.Age5.Witness.KillIncluded = Data.ThreatDepVars$IH5h1 + Data.ThreatDepVars$IH5h3 + Data.ThreatDepVars$IH5h5 + Data.ThreatDepVars$IH5h7
Data.ThreatDepVars$CV.Age5.Witness = Data.ThreatDepVars$IH5h1 + Data.ThreatDepVars$IH5h3 + Data.ThreatDepVars$IH5h5 
Data.ThreatDepVars$CV.Age9.Witness = Data.ThreatDepVars$p5m5a + Data.ThreatDepVars$p5m5b + Data.ThreatDepVars$p5m5c
Data.ThreatDepVars$CV.Age3.Victim = Data.ThreatDepVars$IH3l2 + Data.ThreatDepVars$IH3l4 + Data.ThreatDepVars$IH3l6
Data.ThreatDepVars$CV.Age5.Victim = Data.ThreatDepVars$IH5h2 + Data.ThreatDepVars$IH5h4 + Data.ThreatDepVars$IH5h6
Data.ThreatDepVars$Age3.CV = Data.ThreatDepVars$CV.Age3.Witness.KillIncluded + Data.ThreatDepVars$CV.Age3.Victim
Data.ThreatDepVars$Age5.CV = Data.ThreatDepVars$CV.Age5.Witness.KillIncluded + Data.ThreatDepVars$CV.Age5.Victim
Data.ThreatDepVars$Age9.CV = Data.ThreatDepVars$CV.Age9.Witness

##Neglect ###Item selection based on Hunt, Slack, & Berger 2016 ###Convert home visit neglect items with 7 (yes but not in last year) and 8 (this has never happened) to 0

Data.ThreatDepVars$IH3j15<-car::recode(Data.ThreatDepVars$IH3j15 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH3j16<-car::recode(Data.ThreatDepVars$IH3j16 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH3j17<-car::recode(Data.ThreatDepVars$IH3j17 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH3j18<-car::recode(Data.ThreatDepVars$IH3j18 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH3j19<-car::recode(Data.ThreatDepVars$IH3j19 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH5g15<-car::recode(Data.ThreatDepVars$IH5g15 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH5g16<-car::recode(Data.ThreatDepVars$IH5g16 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH5g17<-car::recode(Data.ThreatDepVars$IH5g17 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH5g18<-car::recode(Data.ThreatDepVars$IH5g18 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$IH5g19<-car::recode(Data.ThreatDepVars$IH5g19 , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$p5q2a<-car::recode(Data.ThreatDepVars$p5q2a , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$p5q2b<-car::recode(Data.ThreatDepVars$p5q2b , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$p5q2c<-car::recode(Data.ThreatDepVars$p5q2c , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$p5q2d<-car::recode(Data.ThreatDepVars$p5q2d , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)
Data.ThreatDepVars$p5q2e<-car::recode(Data.ThreatDepVars$p5q2e , "c(7,8) = 0", as.numeric = TRUE, as.factor = FALSE)

###Calculating physical neglect, emotional neglect, and neglect subscales

Data.ThreatDepVars$Age3.PhysicalNeglect = Data.ThreatDepVars$IH3j15 + Data.ThreatDepVars$IH3j17 + Data.ThreatDepVars$IH3j18 + Data.ThreatDepVars$IH3j19
Data.ThreatDepVars$Age3.EmotionalNeglect =  Data.ThreatDepVars$IH3j16
Data.ThreatDepVars$Age3.Neglect = Data.ThreatDepVars$Age3.PhysicalNeglect + Data.ThreatDepVars$Age3.EmotionalNeglect
Data.ThreatDepVars$Age5.PhysicalNeglect = Data.ThreatDepVars$IH5g15 + Data.ThreatDepVars$IH5g17 + Data.ThreatDepVars$IH5g18 + Data.ThreatDepVars$IH5g19
Data.ThreatDepVars$Age5.EmotionalNeglect = Data.ThreatDepVars$IH5g16
Data.ThreatDepVars$Age5.Neglect = Data.ThreatDepVars$Age5.PhysicalNeglect + Data.ThreatDepVars$Age5.EmotionalNeglect
Data.ThreatDepVars$Age9.PhysicalNeglect = Data.ThreatDepVars$p5q2a + Data.ThreatDepVars$p5q2c + Data.ThreatDepVars$p5q2d + Data.ThreatDepVars$p5q2e
Data.ThreatDepVars$Age9.EmotionalNeglect = Data.ThreatDepVars$p5q2b
Data.ThreatDepVars$Age9.Neglect = Data.ThreatDepVars$Age9.PhysicalNeglect + Data.ThreatDepVars$Age9.EmotionalNeglect

##Neighborhood Collective Efficacy - we’re only focusing on social cohesion ###Item selection based on Donnelly et al., 2016 ###Reverse code items worded in the opposite direction such that higher scores indicate less collective efficacy

Data.ThreatDepVars$IH3k2dr = (6 - Data.ThreatDepVars$IH3k2d)
Data.ThreatDepVars$IH3k2er = (6 - Data.ThreatDepVars$IH3k2e)
Data.ThreatDepVars$m4i0n3r = (5 - Data.ThreatDepVars$m4i0n3)
Data.ThreatDepVars$m4i0n4r = (5 - Data.ThreatDepVars$m4i0n4)
Data.ThreatDepVars$p5m3cr = (5 - Data.ThreatDepVars$p5m3c)
Data.ThreatDepVars$p5m3dr = (5 - Data.ThreatDepVars$p5m3d)

###Calculating social cohesion

Data.ThreatDepVars$Age3.SocialCohesion = Data.ThreatDepVars$IH3k2a + Data.ThreatDepVars$IH3k2b + Data.ThreatDepVars$IH3k2dr + Data.ThreatDepVars$IH3k2er
Data.ThreatDepVars$Age5.SocialCohesion = Data.ThreatDepVars$m4i0n1 + Data.ThreatDepVars$m4i0n2 + Data.ThreatDepVars$m4i0n3r + Data.ThreatDepVars$m4i0n4r
Data.ThreatDepVars$Age9.SocialCohesion = Data.ThreatDepVars$p5m3a + Data.ThreatDepVars$p5m3b + Data.ThreatDepVars$p5m3cr + Data.ThreatDepVars$p5m3dr

##Partner support ###Item selection based on Manuel, Bledsoe, and Bellamy 2012 ###Reverse coded 1 item in opposite direction so that higher score indicates less partner support

Data.ThreatDepVars$m3d7cr = (4 - Data.ThreatDepVars$m3d7c)
Data.ThreatDepVars$m3e23cr = (4 - Data.ThreatDepVars$m3e23c)
Data.ThreatDepVars$m4d7cr = (4 - Data.ThreatDepVars$m4d7c)
Data.ThreatDepVars$m4e23cr = (4 - Data.ThreatDepVars$m4e23c)
Data.ThreatDepVars$m5c6cr = (4 - Data.ThreatDepVars$m5c6c)
Data.ThreatDepVars$m5d20cr = (4 - Data.ThreatDepVars$m5d20c)

###Calculating partner support for baby’s father and current partner

Data.ThreatDepVars$Age3.BF.PartSupp = Data.ThreatDepVars$m3d7a + Data.ThreatDepVars$m3d7b + Data.ThreatDepVars$m3d7cr + Data.ThreatDepVars$m3d7d + Data.ThreatDepVars$m3d7k + Data.ThreatDepVars$m3d7l
Data.ThreatDepVars$Age3.CP.PartSupp = Data.ThreatDepVars$m3e23a + Data.ThreatDepVars$m3e23b + Data.ThreatDepVars$m3e23cr + Data.ThreatDepVars$m3e23d + Data.ThreatDepVars$m3e23k + Data.ThreatDepVars$m3e23l
Data.ThreatDepVars$Age5.BF.PartSupp = Data.ThreatDepVars$m4d7a + Data.ThreatDepVars$m4d7b + Data.ThreatDepVars$m4d7cr + Data.ThreatDepVars$m4d7d + Data.ThreatDepVars$m4d7o + Data.ThreatDepVars$m4d7p
Data.ThreatDepVars$Age5.CP.PartSupp = Data.ThreatDepVars$m4e23a + Data.ThreatDepVars$m4e23b + Data.ThreatDepVars$m4e23cr + Data.ThreatDepVars$m4e23d + Data.ThreatDepVars$m4e23o + Data.ThreatDepVars$m4e23p
Data.ThreatDepVars$Age9.BF.PartSupp = Data.ThreatDepVars$m5c6a + Data.ThreatDepVars$m5c6b + Data.ThreatDepVars$m5c6cr + Data.ThreatDepVars$m5c6d + Data.ThreatDepVars$m5c6o + Data.ThreatDepVars$m5c6p
Data.ThreatDepVars$Age9.CP.PartSupp = Data.ThreatDepVars$m5d20a + Data.ThreatDepVars$m5d20b + Data.ThreatDepVars$m5d20cr + Data.ThreatDepVars$m5d20d + Data.ThreatDepVars$m5d20o + Data.ThreatDepVars$m5d20p

###Merge baby’s father, current partner columns

Data.ThreatDepVars$Age3.PartSupp = Data.ThreatDepVars$Age3.BF.PartSupp
Data.ThreatDepVars$Age3.PartSupp[!is.na(Data.ThreatDepVars$Age3.CP.PartSupp)] = Data.ThreatDepVars$Age3.CP.PartSupp[!is.na(Data.ThreatDepVars$Age3.CP.PartSupp)]
Data.ThreatDepVars$Age5.PartSupp = Data.ThreatDepVars$Age5.BF.PartSupp
Data.ThreatDepVars$Age5.PartSupp[!is.na(Data.ThreatDepVars$Age5.CP.PartSupp)] = Data.ThreatDepVars$Age5.CP.PartSupp[!is.na(Data.ThreatDepVars$Age5.CP.PartSupp)]
Data.ThreatDepVars$Age9.PartSupp = Data.ThreatDepVars$Age9.BF.PartSupp
Data.ThreatDepVars$Age9.PartSupp[!is.na(Data.ThreatDepVars$Age9.CP.PartSupp)] = Data.ThreatDepVars$Age9.CP.PartSupp[!is.na(Data.ThreatDepVars$Age9.CP.PartSupp)]

##Recoding IPV data for participants who did not live with mom at least half the time as NA after email with Chris on 11.09.17

Data.ThreatDepVars$IPV.3.Total <- ifelse(Data.ThreatDepVars$age3cust == 0, NA, Data.ThreatDepVars$IPV.3.Total)
Data.ThreatDepVars$IPV.5.Total <- ifelse(Data.ThreatDepVars$age5cust == 0, NA, Data.ThreatDepVars$IPV.5.Total)
Data.ThreatDepVars$IPV.9.Total <- ifelse(Data.ThreatDepVars$age9cust == 0, NA, Data.ThreatDepVars$IPV.9.Total)

##Recoding partner support for participants who did not live with mom at least half the time as NA

Data.ThreatDepVars$Age3.PartSupp <- ifelse(Data.ThreatDepVars$age3cust == 0, NA, Data.ThreatDepVars$Age3.PartSupp)
Data.ThreatDepVars$Age5.PartSupp <- ifelse(Data.ThreatDepVars$age5cust == 0, NA, Data.ThreatDepVars$Age5.PartSupp)
Data.ThreatDepVars$Age9.PartSupp <- ifelse(Data.ThreatDepVars$age9cust == 0, NA, Data.ThreatDepVars$Age9.PartSupp)

##Creating composite measures of threat and deprivation across ages 3, 5, and 9 ###Calculating Z scores

Data.ThreatDepVars$Age3.AbuseZ <- scale(Data.ThreatDepVars$Age3.Abuse, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age5.AbuseZ <- scale(Data.ThreatDepVars$Age5.Abuse, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age9.AbuseZ <- scale(Data.ThreatDepVars$Age9.Abuse, center = TRUE, scale = TRUE)
Data.ThreatDepVars$IPV.3.TotalZ <- scale(Data.ThreatDepVars$IPV.3.Total, center = TRUE, scale = TRUE) 
Data.ThreatDepVars$IPV.5.TotalZ <- scale(Data.ThreatDepVars$IPV.5.Total, center = TRUE, scale = TRUE) 
Data.ThreatDepVars$IPV.9.TotalZ <- scale(Data.ThreatDepVars$IPV.9.Total, center = TRUE, scale = TRUE) 
Data.ThreatDepVars$Age3.CVZ <- scale(Data.ThreatDepVars$Age3.CV, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age5.CVZ <- scale(Data.ThreatDepVars$Age5.CV, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age9.CVZ <- scale(Data.ThreatDepVars$Age9.CV, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age3.NeglectZ <- scale(Data.ThreatDepVars$Age3.Neglect, center = TRUE, scale = TRUE) 
Data.ThreatDepVars$Age5.NeglectZ <- scale(Data.ThreatDepVars$Age5.Neglect, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age9.NeglectZ <- scale(Data.ThreatDepVars$Age9.Neglect, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age3.SocialCohesionZ <- scale(Data.ThreatDepVars$Age3.SocialCohesion, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age5.SocialCohesionZ <- scale(Data.ThreatDepVars$Age5.SocialCohesion, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age9.SocialCohesionZ <- scale(Data.ThreatDepVars$Age9.SocialCohesion, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age3.PartSuppZ <- scale(Data.ThreatDepVars$Age3.PartSupp, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age5.PartSuppZ <- scale(Data.ThreatDepVars$Age5.PartSupp, center = TRUE, scale = TRUE)
Data.ThreatDepVars$Age9.PartSuppZ <- scale(Data.ThreatDepVars$Age9.PartSupp, center = TRUE, scale = TRUE)

###Creating an indicator of how many exposures each participant has data for

Data.ThreatDepVars$Age3.Abuse.Incl = ifelse(is.na(Data.ThreatDepVars$Age3.Abuse), "0", "1")
Data.ThreatDepVars$Age3.Abuse.Incl <- as.numeric(Data.ThreatDepVars$Age3.Abuse.Incl)
Data.ThreatDepVars$Age5.Abuse.Incl = ifelse(is.na(Data.ThreatDepVars$Age5.Abuse), "0", "1")
Data.ThreatDepVars$Age5.Abuse.Incl <- as.numeric(Data.ThreatDepVars$Age5.Abuse.Incl)
Data.ThreatDepVars$Age9.Abuse.Incl = ifelse(is.na(Data.ThreatDepVars$Age9.Abuse), "0", "1")
Data.ThreatDepVars$Age9.Abuse.Incl <- as.numeric(Data.ThreatDepVars$Age9.Abuse.Incl)
Data.ThreatDepVars$Age3.IPV.Incl = ifelse(is.na(Data.ThreatDepVars$IPV.3.Total), "0", "1")
Data.ThreatDepVars$Age3.IPV.Incl <- as.numeric(Data.ThreatDepVars$Age3.IPV.Incl)
Data.ThreatDepVars$Age5.IPV.Incl = ifelse(is.na(Data.ThreatDepVars$IPV.5.Total), "0", "1")
Data.ThreatDepVars$Age5.IPV.Incl <- as.numeric(Data.ThreatDepVars$Age5.IPV.Incl)
Data.ThreatDepVars$Age9.IPV.Incl = ifelse(is.na(Data.ThreatDepVars$IPV.9.Total), "0", "1")
Data.ThreatDepVars$Age9.IPV.Incl <- as.numeric(Data.ThreatDepVars$Age9.IPV.Incl)
Data.ThreatDepVars$Age3.CV.Incl = ifelse(is.na(Data.ThreatDepVars$Age3.CV), "0", "1")
Data.ThreatDepVars$Age3.CV.Incl <- as.numeric(Data.ThreatDepVars$Age3.CV.Incl)
Data.ThreatDepVars$Age5.CV.Incl = ifelse(is.na(Data.ThreatDepVars$Age5.CV), "0", "1")
Data.ThreatDepVars$Age5.CV.Incl <- as.numeric(Data.ThreatDepVars$Age5.CV.Incl)
Data.ThreatDepVars$Age9.CV.Incl = ifelse(is.na(Data.ThreatDepVars$Age9.CV), "0", "1")
Data.ThreatDepVars$Age9.CV.Incl <- as.numeric(Data.ThreatDepVars$Age9.CV.Incl)
Data.ThreatDepVars$Age3.Neglect.Incl = ifelse(is.na(Data.ThreatDepVars$Age3.Neglect), "0", "1")
Data.ThreatDepVars$Age3.Neglect.Incl <- as.numeric(Data.ThreatDepVars$Age3.Neglect.Incl)
Data.ThreatDepVars$Age5.Neglect.Incl = ifelse(is.na(Data.ThreatDepVars$Age5.Neglect), "0", "1")
Data.ThreatDepVars$Age5.Neglect.Incl <- as.numeric(Data.ThreatDepVars$Age5.Neglect.Incl)
Data.ThreatDepVars$Age9.Neglect.Incl = ifelse(is.na(Data.ThreatDepVars$Age9.Neglect), "0", "1")
Data.ThreatDepVars$Age9.Neglect.Incl <- as.numeric(Data.ThreatDepVars$Age9.Neglect.Incl)
Data.ThreatDepVars$Age3.SocialCohesion.Incl = ifelse(is.na(Data.ThreatDepVars$Age3.SocialCohesion), "0", "1")
Data.ThreatDepVars$Age3.SocialCohesion.Incl <- as.numeric(Data.ThreatDepVars$Age3.SocialCohesion.Incl)
Data.ThreatDepVars$Age5.SocialCohesion.Incl = ifelse(is.na(Data.ThreatDepVars$Age5.SocialCohesion), "0", "1")
Data.ThreatDepVars$Age5.SocialCohesion.Incl <- as.numeric(Data.ThreatDepVars$Age5.SocialCohesion.Incl)
Data.ThreatDepVars$Age9.SocialCohesion.Incl = ifelse(is.na(Data.ThreatDepVars$Age9.SocialCohesion), "0", "1")
Data.ThreatDepVars$Age9.SocialCohesion.Incl <- as.numeric(Data.ThreatDepVars$Age9.SocialCohesion.Incl)
Data.ThreatDepVars$Age3.PartSupp.Incl = ifelse(is.na(Data.ThreatDepVars$Age3.PartSupp), "0", "1")
Data.ThreatDepVars$Age3.PartSupp.Incl <- as.numeric(Data.ThreatDepVars$Age3.PartSupp.Incl)
Data.ThreatDepVars$Age5.PartSupp.Incl = ifelse(is.na(Data.ThreatDepVars$Age5.PartSupp), "0", "1")
Data.ThreatDepVars$Age5.PartSupp.Incl <- as.numeric(Data.ThreatDepVars$Age5.PartSupp.Incl)
Data.ThreatDepVars$Age9.PartSupp.Incl = ifelse(is.na(Data.ThreatDepVars$Age9.PartSupp), "0", "1")
Data.ThreatDepVars$Age9.PartSupp.Incl <- as.numeric(Data.ThreatDepVars$Age9.PartSupp.Incl)

Data.ThreatDepVars$Threat.Incl = Data.ThreatDepVars$Age3.Abuse.Incl + Data.ThreatDepVars$Age5.Abuse.Incl + Data.ThreatDepVars$Age9.Abuse.Incl + Data.ThreatDepVars$Age3.IPV.Incl + Data.ThreatDepVars$Age5.IPV.Incl + Data.ThreatDepVars$Age9.IPV.Incl + Data.ThreatDepVars$Age3.CV.Incl + Data.ThreatDepVars$Age5.CV.Incl +  Data.ThreatDepVars$Age9.CV.Incl
Data.ThreatDepVars$Dep.Incl = Data.ThreatDepVars$Age3.Neglect.Incl + Data.ThreatDepVars$Age5.Neglect.Incl + Data.ThreatDepVars$Age9.Neglect.Incl + Data.ThreatDepVars$Age3.SocialCohesion.Incl + Data.ThreatDepVars$Age5.SocialCohesion.Incl + Data.ThreatDepVars$Age9.SocialCohesion.Incl + Data.ThreatDepVars$Age3.PartSupp.Incl + Data.ThreatDepVars$Age5.PartSupp.Incl + Data.ThreatDepVars$Age9.PartSupp.Incl

###Creating composite scores for threat and deprivation

Data.ThreatDepVars$ThreatZSum <- rowSums(cbind(Data.ThreatDepVars$Age3.AbuseZ,Data.ThreatDepVars$Age5.AbuseZ,Data.ThreatDepVars$Age9.AbuseZ,Data.ThreatDepVars$IPV.3.TotalZ,Data.ThreatDepVars$IPV.5.TotalZ,Data.ThreatDepVars$IPV.9.TotalZ,Data.ThreatDepVars$Age3.CVZ,Data.ThreatDepVars$Age5.CVZ,Data.ThreatDepVars$Age9.CVZ), na.rm = TRUE)
Data.ThreatDepVars$ThreatComp <- Data.ThreatDepVars$ThreatZSum/Data.ThreatDepVars$Threat.Incl
Data.ThreatDepVars$DepZSum <-rowSums(cbind(Data.ThreatDepVars$Age3.NeglectZ,Data.ThreatDepVars$Age5.NeglectZ,Data.ThreatDepVars$Age9.NeglectZ,Data.ThreatDepVars$Age3.SocialCohesionZ,Data.ThreatDepVars$Age5.SocialCohesionZ,Data.ThreatDepVars$Age9.SocialCohesionZ,Data.ThreatDepVars$Age3.PartSuppZ,Data.ThreatDepVars$Age5.PartSuppZ,Data.ThreatDepVars$Age9.PartSuppZ), na.rm=TRUE)
Data.ThreatDepVars$DepComp <- Data.ThreatDepVars$DepZSum/Data.ThreatDepVars$Dep.Incl

Select only the Composites

Data.ThreatDepVars = Data.ThreatDepVars %>% 
  dplyr::select("ff_id", "ThreatComp", "DepComp")

Demographic Variable Prep

# Race/Ethnicity at Age 15- ck6ethrace
# 5 Multi-racial, non-hispanic
# 4 Other only, non-hispanic
# 3 Hispanic/Latino
# 2 Black/Af. American only, non-hispanic
# 1 White only, non-hispanic

Demo_Variables = Demo_Variables %>% 
  rowwise() %>% 
  mutate(povco_sum = sum(cm1inpov, cm2povco, cm3povco, cm4povco, cm5povco, cp6povco, na.rm=TRUE)) %>% 
  mutate(wv1_pov_inc = if_else(is.na(cm1inpov), 0,1),
         wv2_pov_inc = if_else(is.na(cm2povco), 0,1),
         wv3_pov_inc = if_else(is.na(cm3povco), 0,1),
         wv4_pov_inc = if_else(is.na(cm4povco), 0,1),
         wv5_pov_inc = if_else(is.na(cm5povco), 0,1),
         wv6_pov_inc = if_else(is.na(cp6povco), 0,1)) %>% 
  mutate(pov_inc_sum = sum(wv1_pov_inc,wv2_pov_inc,wv3_pov_inc,wv4_pov_inc,wv5_pov_inc,wv6_pov_inc, na.rm = TRUE)) %>% 
  mutate(povco_avg = povco_sum/pov_inc_sum) %>% 
  mutate(Race_AA = if_else(ck6ethrace == 2, 1, 0),
         Race_C = if_else(ck6ethrace == 1, 1, 0),
         Race_L = if_else(ck6ethrace == 3, 1, 0)) %>% 
  mutate(cm1bsex = cm1bsex-1) %>% 
  dplyr::select(ff_id, povco_avg, Race_AA, Race_C, Race_L, ck6ethrace, cm1bsex, m1city)

Merge to One Dataframe

All_Variables = Data.ThreatDepVars %>% 
  left_join(Data.AnxDep, by="ff_id")

All_Variables = All_Variables %>% 
  left_join(Extern_Variables, by="ff_id")

All_Variables = All_Variables %>% 
  left_join(Data.PAF, by = "ff_id")

All_Variables = All_Variables %>% 
  left_join(Data.SchoolConn, by = "ff_id")

All_Variables = All_Variables %>% 
  left_join(CBCL, by="ff_id")

All_Variables = All_Variables %>% 
  left_join(Demo_Variables, by = "ff_id")

Write Variables to CSV

write.csv(All_Variables,'/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_012021.csv')

Write Variables for MPlus

All_Variables[is.na(All_Variables)] = 99
prepareMplusData(All_Variables,'All_Variables_012021.dat')

Recode Missing Back to NA for Further Use

All_Variables[is.na(All_Variables)] = 99

Actual Analysis Starts Here - Read in Data

Here I am reading in data written out from the creation/recoding of variables above as well as the factor scores extracted from Mplus for the simple slopes analysis. Convention for how factor scores should be extracted when probing latent variable interactions has not been clearly established. Therefore, we used two different approaches and tested whether the method used influenced the interpretation of the interaction. The first method extracted the factor scores from a measurement model containing all of the latent variables in the model (e.g. school connectedness at age 9, school connectedness at age 15, and positive adolescent function). The second approach used the factor scores extracted from individual measurement models. Due to the correlations between latent variables that are accounted for in measurement models (Kline, 2015), it is likely that the extracted factors are different depending on how they are extracted. This is why there are two dataframes, All_Variables and All_Variables2.

All_Variables = read_csv('VE_SD_IntExtPAF_SchoolConn_012021.csv')
## Parsed with column specification:
## cols(
##   .default = col_double()
## )
## See spec(...) for full column specifications.
# List of IDs that are included in models with covariates
fs_IDs = read_csv('ff_ids.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   city = col_double()
## )
fs_IDs$ff_id = as.integer(fs_IDs$ff_id)

# Factor scores from SC9,SC15,PAF Measurement Mod (stratified)
fs_scores_pafmod = read_csv('factor_scores_SC15_SC9_PAF_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_character(),
##   SC15 = col_double(),
##   SC15_SE = col_double(),
##   SC9 = col_double(),
##   SC9_SE = col_double(),
##   PAF = col_double(),
##   PAF_SE = col_double(),
##   city = col_double()
## )
fs_scores_pafmod$ff_id = as.integer(fs_scores_pafmod$ff_id)
## Warning: NAs introduced by coercion
# Factor scores from SC9,Ext9,Int9 Measurement Mod (stratified)
fs_scores_ext9mod = read_csv('factor_scores_SC_Int_Ext_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_character(),
##   SC15 = col_double(),
##   SC15_SE = col_double(),
##   SC9 = col_double(),
##   SC9_SE = col_double(),
##   Int15 = col_double(),
##   Int15_SE = col_double(),
##   Ext15 = col_double(),
##   Ext15_SE = col_double(),
##   Int9 = col_double(),
##   Int9_SE = col_double(),
##   Ext9 = col_double(),
##   Ext9_SE = col_double(),
##   city = col_double()
## )
fs_scores_ext9mod$ff_id = as.integer(fs_scores_ext9mod$ff_id)
## Warning: NAs introduced by coercion
# Individual factor scores -- needs to be updated to stratified
SC9 = read_csv('factor_scores_SC9_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_character(),
##   SC9 = col_double(),
##   SC9_SE = col_double(),
##   city = col_double()
## )
SC9$ff_id = as.integer(SC9$ff_id)
## Warning: NAs introduced by coercion
SC15 = read_csv('factor_scores_SC15_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   SC15 = col_double(),
##   SC15_SE = col_double(),
##   city = col_double()
## )
SC15$ff_id = as.integer(SC15$ff_id)

Int9 = read_csv('factor_scores_Int9_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_character(),
##   INT9 = col_double(),
##   INT9_SE = col_double(),
##   city = col_double()
## )
Int9$ff_id = as.integer(Int9$ff_id)
## Warning: NAs introduced by coercion
Ext9 = read_csv('factor_scores_Ext9_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_character(),
##   EXT9 = col_double(),
##   EXT9_SE = col_double(),
##   city = col_double()
## )
Ext9$ff_id = as.integer(Ext9$ff_id)
## Warning: NAs introduced by coercion
Int15 = read_csv('factor_scores_Int15_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   INT15 = col_double(),
##   INT15_SE = col_double(),
##   city = col_double()
## )
Int15$ff_id = as.integer(Int15$ff_id)

Ext15 = read_csv('factor_scores_Ext15_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   EXT15 = col_double(),
##   EXT15_SE = col_double(),
##   city = col_double()
## )
Ext15$ff_id = as.integer(Ext15$ff_id)

PAF = read_csv('factor_scores_PAF_012221.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   PAF = col_double(),
##   PAF_SE = col_double(),
##   city = col_double()
## )
PAF$ff_id = as.integer(PAF$ff_id)

# Merge things

PAF_mod_data = fs_IDs %>% 
  left_join(fs_scores_pafmod) %>% 
  left_join(All_Variables) 
## Joining, by = c("ff_id", "city")
## Joining, by = "ff_id"
Ext9_mod_data = fs_IDs %>% 
  left_join(fs_scores_ext9mod) %>% 
  left_join(All_Variables) 
## Joining, by = c("ff_id", "city")
## Joining, by = "ff_id"
Indiv_FS_All_Variables = fs_IDs %>%
  left_join(All_Variables) %>% 
  left_join(SC9) %>% 
  left_join(SC15) %>% 
  left_join(Int9) %>% 
  left_join(Ext9) %>% 
  left_join(Int15) %>% 
  left_join(Ext15) %>% 
  left_join(PAF)
## Joining, by = "ff_id"
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")

This reverses any 99s that are marked in the MPlus data – this should already be done, but just in case.

is.na(PAF_mod_data) = PAF_mod_data == 99
is.na(Ext9_mod_data) = Ext9_mod_data == 99
is.na(Indiv_FS_All_Variables) = Indiv_FS_All_Variables == 99

Look at distributions of our variables

PAF = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=PAF)) + xlab("Positive Adolescent Function") + geom_density() + theme_classic()

Dep = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=DepComp)) + xlab("Social Deprivation") + geom_density() + theme_classic()

Threat = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=ThreatComp)) + xlab("Violence Exposure") + geom_density() + theme_classic()

SC9 = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=SC9)) + xlab("School Connectedness Age 9") + geom_density() + theme_classic()

SC15 = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=SC15)) + xlab("School Connectedness Age 15") + geom_density() + theme_classic()

Int9 =  Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=INT9)) + xlab("Internalizing Psychopathology Age 9") + geom_density() + theme_classic()

Ext9 =  Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=EXT9)) + xlab("Externalizing Psychopathology Age 9") + geom_density() + theme_classic()

Int15 =  Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=INT15)) + xlab("Internalizing Psychopathology Age 15") + geom_density() + theme_classic()

Ext15 = Indiv_FS_All_Variables %>% 
  ggplot(mapping = aes(x=EXT15)) + xlab("Externalizing Psychopathology Age 15") + geom_density() + theme_classic()

figure <- ggarrange(PAF, Threat, Dep, SC9,SC15, Int9, Ext9, Int15,Ext15,
                    ncol = 2, nrow = 5)
## Warning: Removed 2 rows containing non-finite values (stat_density).
## Warning: Removed 353 rows containing non-finite values (stat_density).
## Warning: Removed 53 rows containing non-finite values (stat_density).
## Warning: Removed 366 rows containing non-finite values (stat_density).

## Warning: Removed 366 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing non-finite values (stat_density).
## Warning: Removed 1 rows containing non-finite values (stat_density).
figure

Skewness and Kurtosis

psych::describe(Indiv_FS_All_Variables$ThreatComp) # -- this is potentially problematic, it seems like there are some serious kurtosis issues 
##    vars    n mean   sd median trimmed  mad   min max range skew kurtosis   se
## X1    1 3246 0.01 0.53  -0.09   -0.05 0.45 -1.13 7.1  8.23 1.95    11.84 0.01
psych::describe(Indiv_FS_All_Variables$DepComp)
##    vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 3246    0 0.53  -0.08   -0.04 0.47 -1.47 4.02  5.49 1.34      4.3 0.01
psych::describe(Indiv_FS_All_Variables$PAF)
##    vars    n  mean   sd median trimmed  mad   min  max range  skew kurtosis
## X1    1 3244 -0.02 0.89  -0.01   -0.01 0.93 -3.04 1.56   4.6 -0.15    -0.39
##      se
## X1 0.02
psych::describe(Indiv_FS_All_Variables$INT9)
##    vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 2880 0.08 0.82   0.03    0.02 0.86 -1.01 4.92  5.92 0.73     1.38 0.02
psych::describe(Indiv_FS_All_Variables$EXT9)
##    vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 2880 0.05 0.88   0.02       0 0.95 -1.21 4.28  5.48 0.44     0.06 0.02
psych::describe(Indiv_FS_All_Variables$INT15)
##    vars    n mean   sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 3244 0.03 0.87   0.02       0 0.95 -1.37 3.18  4.55 0.24    -0.46 0.02
psych::describe(Indiv_FS_All_Variables$EXT15)
##    vars    n mean   sd median trimmed mad   min max range  skew kurtosis   se
## X1    1 3245 0.04 0.91   0.07    0.06 0.9 -2.06 2.8  4.86 -0.13    -0.22 0.02
psych::describe(Indiv_FS_All_Variables$SC9)
##    vars    n  mean   sd median trimmed mad   min  max range  skew kurtosis   se
## X1    1 2893 -0.05 0.77   0.01    0.02 0.9 -2.28 0.86  3.14 -0.59     -0.2 0.01
psych::describe(Indiv_FS_All_Variables$SC15)
##    vars    n  mean   sd median trimmed  mad   min  max range  skew kurtosis
## X1    1 3193 -0.04 0.78  -0.02    0.03 0.76 -2.65 0.91  3.57 -0.56    -0.22
##      se
## X1 0.01

Table Creation

The code below creates: - Demographic Table - Correlation Table w/descriptives for our main variables.

# Import a categorical race variable because I had accidentally cut it out earlier.

Race_Cat = read_csv('Race_Categorical.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   Race_Cat = col_character()
## )
Race_Cat$Race_Cat = as.factor(Race_Cat$Race_Cat)

MatBirthVars = read_csv('Maternal_Birth_Variables.csv')
## Parsed with column specification:
## cols(
##   ff_id = col_double(),
##   m1b2 = col_double(),
##   cm1edu = col_double()
## )
MatBirthVars$m1b2 = as.factor(MatBirthVars$m1b2)
MatBirthVars$cm1edu = as.factor(MatBirthVars$cm1edu)

Indiv_FS_All_Variables = Indiv_FS_All_Variables %>% 
  left_join(Race_Cat, by='ff_id') %>% 
  left_join(MatBirthVars, by='ff_id')


Indiv_FS_All_Variables = Indiv_FS_All_Variables %>% 
  mutate(MatMar = case_when(
    m1b2 == 1 ~ 'Married',
    m1b2 == 2 ~ 'Not Married',
    TRUE ~ 'NA'
  ),
  MatEdu = case_when(
    cm1edu == 1 ~ 'Less than high school',
    cm1edu == 2 ~ 'High school or equivalent',
    cm1edu == 3 ~ 'Some college or technical school',
    cm1edu == 4 ~ 'College or graduate school',
    TRUE ~ 'NA'
  ),
  sex = case_when(
    cm1bsex == 0 ~ 'Male',
    cm1bsex == 1 ~ 'Female',
    TRUE ~ 'NA'
  ))
table1::label(Indiv_FS_All_Variables$sex) = "Sex"
table1::label(Indiv_FS_All_Variables$ck6ethrace) = "Race-Ethnicity"
table1::label(Indiv_FS_All_Variables$povco_avg) = "Average Poverty Ratio"
table1::label(Indiv_FS_All_Variables$MatMar) = "Maternal Marital Status at Child's Birth"
table1::label(Indiv_FS_All_Variables$MatEdu) = "Maternal Education at Child's Birth"

table1::label(Indiv_FS_All_Variables$ThreatComp) = "Violence Exposure"
table1::label(Indiv_FS_All_Variables$DepComp) = "Social Deprivation"
table1::label(Indiv_FS_All_Variables$SC9) = "School Connectedness Age 9"
table1::label(Indiv_FS_All_Variables$SC15) = "School Connectedness Age 15"
table1::label(Indiv_FS_All_Variables$PAF) = "Positive Adolescent Function"
table1::label(Indiv_FS_All_Variables$INT15) = "Internalizing Symptoms Age 15"
table1::label(Indiv_FS_All_Variables$EXT15) = "Externalizing Symptoms Age 15"
table1::label(Indiv_FS_All_Variables$INT9) = "Internalizing Symptoms Age 9"
table1::label(Indiv_FS_All_Variables$EXT9) = "Externalizing Symptoms Age 9"

table1::table1(~sex + Race_Cat + povco_avg + MatMar + MatEdu, data = Indiv_FS_All_Variables)
Overall
(N=3246)
Sex
Female 1585 (48.8%)
Male 1661 (51.2%)
Race_Cat
African American 1592 (49.0%)
Caucasian 587 (18.1%)
Latinx 808 (24.9%)
Other 259 (8.0%)
Average Poverty Ratio
Mean (SD) 2.11 (2.10)
Median [Min, Max] 1.46 [0.120, 21.2]
Maternal Marital Status at Child's Birth
Married 785 (24.2%)
NA 18 (0.6%)
Not Married 2443 (75.3%)
Maternal Education at Child's Birth
College or graduate school 365 (11.2%)
High school or equivalent 1030 (31.7%)
Less than high school 1025 (31.6%)
NA 5 (0.2%)
Some college or technical school 821 (25.3%)
table1::table1(~ThreatComp + DepComp + SC9 + SC15 + INT9 + EXT9 + INT15 + EXT15 + PAF, data = Indiv_FS_All_Variables)
Overall
(N=3246)
Violence Exposure
Mean (SD) 0.00869 (0.535)
Median [Min, Max] -0.0922 [-1.13, 7.10]
Social Deprivation
Mean (SD) 0.00323 (0.533)
Median [Min, Max] -0.0780 [-1.47, 4.02]
School Connectedness Age 9
Mean (SD) -0.0524 (0.769)
Median [Min, Max] 0.0100 [-2.28, 0.865]
Missing 353 (10.9%)
School Connectedness Age 15
Mean (SD) -0.0440 (0.784)
Median [Min, Max] -0.0190 [-2.66, 0.912]
Missing 53 (1.6%)
Internalizing Symptoms Age 9
Mean (SD) 0.0752 (0.821)
Median [Min, Max] 0.0310 [-1.01, 4.92]
Missing 366 (11.3%)
Externalizing Symptoms Age 9
Mean (SD) 0.0467 (0.878)
Median [Min, Max] 0.0230 [-1.21, 4.28]
Missing 366 (11.3%)
Internalizing Symptoms Age 15
Mean (SD) 0.0257 (0.873)
Median [Min, Max] 0.0185 [-1.37, 3.18]
Missing 2 (0.1%)
Externalizing Symptoms Age 15
Mean (SD) 0.0418 (0.908)
Median [Min, Max] 0.0710 [-2.06, 2.80]
Missing 1 (0.0%)
Positive Adolescent Function
Mean (SD) -0.0182 (0.894)
Median [Min, Max] -0.0140 [-3.04, 1.56]
Missing 2 (0.1%)
CorTable = Indiv_FS_All_Variables %>% 
  dplyr::select(ThreatComp, DepComp, SC9, SC15, INT9, EXT9, INT15, EXT15, PAF)

apa.cor.table(CorTable, filename = '/Users/leighgoetschius/Box Sync/Shift_Persist_FF/Manuscript/Figures/cor_table_012321.doc', table.number = NA,
  show.conf.interval = TRUE, landscape = TRUE)

Quick check of VE-SD associations

cor.test(Indiv_FS_All_Variables$ThreatComp, Indiv_FS_All_Variables$DepComp)
## 
##  Pearson's product-moment correlation
## 
## data:  Indiv_FS_All_Variables$ThreatComp and Indiv_FS_All_Variables$DepComp
## t = 24.925, df = 3244, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3716339 0.4293926
## sample estimates:
##       cor 
## 0.4009115
fit_cor = lm(PAF ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor2 = lm(EXT9 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor3 = lm(INT9 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor4 = lm(EXT15 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor5 = lm(INT15 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)

paf_vif = car::vif(fit_cor)
ext9_vif = car::vif(fit_cor2)
int9_vif = car::vif(fit_cor3)
ext15_vif = car::vif(fit_cor4)
int15_vif = car::vif(fit_cor5)

mean(paf_vif[1],ext9_vif[1],int9_vif[1],ext15_vif[1],int15_vif[1])
## [1] 1.190807

Plot of the SC15 Effect

PAF_form = lm(Indiv_FS_All_Variables$PAF~Indiv_FS_All_Variables$SC15)
Int_form = lm(Indiv_FS_All_Variables$INT15~Indiv_FS_All_Variables$SC15)
Ext_form = lm(Indiv_FS_All_Variables$EXT15~Indiv_FS_All_Variables$SC15)

PAF_SC15 = Indiv_FS_All_Variables %>% 
  ggplot(aes(x=SC15, y=PAF)) +
  geom_point(position="jitter") +
  geom_smooth(method="lm", colour="skyblue3") + 
  xlab("School Connectedness Age 15") + 
  ylab("Positive Adolescent Function Age 15") +
  stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = PAF_form) +
  theme_classic() + 
  theme(text = element_text(family='serif'))

Int_SC15 = Indiv_FS_All_Variables %>% 
  ggplot(aes(x=SC15, y=INT15)) +
  geom_point(position="jitter") +
  geom_smooth(method="lm", colour="red") +
  xlab("School Connectedness Age 15") + 
  ylab("Internalizing Symptoms Age 15") +
  stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int_form) +
  theme_classic() + 
  theme(text = element_text(family='serif'))


Ext_SC15 = Indiv_FS_All_Variables %>% 
  ggplot(aes(x=SC15, y=EXT15)) +
  geom_point(position="jitter") + 
  geom_smooth(method="lm", colour = "red")+
  xlab("School Connectedness Age 15") + 
  ylab("Externalizing Symptoms Age 15") +
  stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext_form) +
  theme_classic() + 
  theme(text = element_text(family='serif'))



SC15_figure <- ggarrange(Int_SC15, Ext_SC15, PAF_SC15,
                    ncol = 2, nrow = 2)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 55 rows containing non-finite values (stat_smooth).
## Warning: Removed 55 rows containing non-finite values (stat_regline_equation).
## Warning: Removed 55 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 54 rows containing non-finite values (stat_smooth).
## Warning: Removed 54 rows containing non-finite values (stat_regline_equation).
## Warning: Removed 54 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 55 rows containing non-finite values (stat_smooth).
## Warning: Removed 55 rows containing non-finite values (stat_regline_equation).
## Warning: Removed 55 rows containing missing values (geom_point).
SC15_figure

Plot SC9 Contemporaneous Effect

Int9_form = lm(Indiv_FS_All_Variables$INT9~Indiv_FS_All_Variables$SC9)
Ext9_form = lm(Indiv_FS_All_Variables$EXT9~Indiv_FS_All_Variables$SC9)

Int9_SC9 = Indiv_FS_All_Variables %>% 
  ggplot(aes(x=SC9, y=INT9)) +
  geom_point(position="jitter") +
  geom_smooth(method="lm", colour="red") +
  xlab("School Connectedness Age 9") + 
  ylab("Internalizing Symptoms Age 9") +
  stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int9_form) +
  theme_classic() + 
  theme(text = element_text(size=20, family='serif'))


Ext9_SC9 = Indiv_FS_All_Variables %>% 
  ggplot(aes(x=SC9, y=EXT9)) +
  geom_point(position="jitter") + 
  geom_smooth(method="lm", colour = "red")+
  xlab("School Connectedness Age 9") + 
  ylab("Externalizing Symptoms Age 9") +
  theme_classic() + 
  theme(text = element_text(size=20, family='serif')) +
  stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext9_form) 



SC9_figure <- ggarrange(Int9_SC9, Ext9_SC9,
                    ncol = 2, nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 399 rows containing non-finite values (stat_smooth).
## Warning: Removed 399 rows containing non-finite values (stat_regline_equation).
## Warning: Removed 399 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 399 rows containing non-finite values (stat_smooth).
## Warning: Removed 399 rows containing non-finite values (stat_regline_equation).
## Warning: Removed 399 rows containing missing values (geom_point).
SC9_figure

Probing SC9xSD –> PAF w.stratification

This uses factor scores extracted from a measurement model with SC9, SC15, and PAF that is stratified by m1city (city at birth)

SC9 as moderator

# Center/Standardize the Variables
PAF_mod_data = PAF_mod_data %>% 
  mutate(PAF_scaled = scale(PAF),
         DepComp_scaled = scale(DepComp),
         SC9_scaled = scale(SC9))

# Set your model to a variable
sc9xpaf_mod = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data)

sc9PAF_ros = data.frame(xmin = -2.772,
                       xmax = 2.766,
                       ymin = -Inf, 
                       ymax = Inf)

# Interpret the interaction
interact_plot(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "Social Deprivation", y.label = "Positive Function Age 15", legend.main = "School Connectedness", data = PAF_mod_data) + geom_rect(data=sc9PAF_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.3, inherit.aes=FALSE)  + geom_vline(xintercept = 0, colour="black", linetype = "longdash") + theme_classic() + theme(text=element_text(size = 16, family="serif"))

johnson_neyman(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is OUTSIDE the interval [-15.914, -1.270], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-3.262, 1.456]

johnson_neyman(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is OUTSIDE the interval [2.766, 32.588], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]

sim_slopes(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, data = PAF_mod_data, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is OUTSIDE the interval [-15.914, -1.270], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-3.262, 1.456]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of DepComp_scaled when SC9_scaled = -1.000 (- 1 SD): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.066   0.024   -2.729   0.006
## 
## Slope of DepComp_scaled when SC9_scaled = -0.000 (Mean): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.109   0.017   -6.381   0.000
## 
## Slope of DepComp_scaled when SC9_scaled =  1.000 (+ 1 SD): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.151   0.025   -5.949   0.000

SD as moderator

SDPAF_ros = data.frame(xmin = -1.270,
                       xmax = 1.456,
                       ymin = -Inf, 
                       ymax = Inf)

interact_plot(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "School Connectedness (Age 9)", y.label = "Positive Function Age 15", legend.main = "Social Deprivation", data = PAF_mod_data) + geom_rect(data=SDPAF_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.2, inherit.aes=FALSE)  + geom_vline(xintercept = 0, colour="black", linetype = "longdash") + theme_classic() + theme(text=element_text(size = 20, family="serif"))

sim_slopes(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, data = PAF_mod_data, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is OUTSIDE the interval [2.766, 32.588], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of SC9_scaled when DepComp_scaled = -1.000 (- 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.263   0.025   10.636   0.000
## 
## Slope of SC9_scaled when DepComp_scaled = -0.000 (Mean): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.221   0.017   12.981   0.000
## 
## Slope of SC9_scaled when DepComp_scaled =  1.000 (+ 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.179   0.025    7.160   0.000

Getting coefficients, variances, and covariances for and running the Preacher interaction tool.

vcov(sc9xpaf_mod)
##                            (Intercept) DepComp_scaled   SC9_scaled
## (Intercept)               2.896838e-04   8.515363e-07 1.622477e-07
## DepComp_scaled            8.515363e-07   2.901365e-04 1.892200e-05
## SC9_scaled                1.622477e-07   1.892200e-05 2.896314e-04
## DepComp_scaled:SC9_scaled 2.129388e-05   1.310645e-05 2.497242e-06
##                           DepComp_scaled:SC9_scaled
## (Intercept)                            2.129388e-05
## DepComp_scaled                         1.310645e-05
## SC9_scaled                             2.497242e-06
## DepComp_scaled:SC9_scaled              3.277457e-04
summary(sc9xpaf_mod)
## 
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3807 -0.6616  0.0105  0.6985  2.2747 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -0.002748   0.017020  -0.161   0.8717    
## DepComp_scaled            -0.108682   0.017033  -6.381 2.02e-10 ***
## SC9_scaled                 0.220912   0.017019  12.981  < 2e-16 ***
## DepComp_scaled:SC9_scaled -0.042298   0.018104  -2.336   0.0195 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9674 on 3242 degrees of freedom
## Multiple R-squared:  0.06504,    Adjusted R-squared:  0.06418 
## F-statistic: 75.18 on 3 and 3242 DF,  p-value: < 2.2e-16
psych::describe(PAF_mod_data$DepComp_scaled)
##    vars    n mean sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 3246    0  1  -0.15   -0.09 0.89 -2.77 7.54 10.31 1.34      4.3 0.02
psych::describe(PAF_mod_data$SC9_scaled)
##    vars    n mean sd median trimmed mad   min  max range  skew kurtosis   se
## X1    1 3246    0  1   0.07    0.08   1 -3.26 1.46  4.72 -0.63     0.11 0.02

Preacher Output - SC9 as Moderator

Your Input

X1 = -2.77 X2 = 7.54 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = -0.002748 X Slope = -0.108682 Z Slope = 0.220912 XZ Slope = -0.042298 df = 3240 alpha = 0.05

Asymptotic (Co)variances

var(b0) 0.00028968 var(b1) 0.00029014 var(b2) 0.00028963 var(b3) 0.00032775 cov(b2,b0) 1.6e-7 cov(b3,b1) 0.00001311

Region of Significance

Z at lower bound of region = -15.914 Z at upper bound of region = -1.2702 (simple slopes are significant outside this region.)

Simple Intercepts and Slopes at Conditional Values of Z

At Z = cv1… simple intercept = -0.2237(0.0241), t=-9.2951, p=0 simple slope = -0.0664(0.0243), t=-2.7291, p=0.0064 At Z = cv2… simple intercept = -0.0027(0.017), t=-0.1615, p=0.8717 simple slope = -0.1087(0.017), t=-6.3805, p=0 At Z = cv3… simple intercept = 0.2182(0.0241), t=9.0616, p=0 simple slope = -0.151(0.0254), t=-5.949, p=0

Simple Intercepts and Slopes at Region Boundaries

Lower Bound…
simple intercept = -3.5183(0.2714), t=-12.9657, p=0 simple slope = 0.5644(0.2879), t=1.9607, p=0.05 Upper Bound…
simple intercept = -0.2834(0.0275), t=-10.3015, p=0 simple slope = -0.055(0.028), t=-1.9606, p=0.05

Points to Plot

Line for cv1: From {X=-2.77, Y=-0.0398} to {X=7.54, Y=-0.7242} Line for cv2: From {X=-2.77, Y=0.2983} to {X=7.54, Y=-0.8222} Line for cv3: From {X=-2.77, Y=0.6364} to {X=7.54, Y=-0.9202}

Preacher Simple Slopes Plot

xx <- c(-2.77,7.54)   #  <-- change to alter plot dims
yy <- c(-1.2315,0.6364)   #  <-- change to alter plot dims
leg <- c(-2.77,-0.9981)   #  <-- change to alter legend location
x <- c(-2.77,7.54)   #  <-- x-coords for lines
y1 <- c(-0.0398,-0.7242)
y2 <- c(0.2983,-0.8222)
y3 <- c(0.6364,-0.9202)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='X',ylab='Y',main='MLR 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

Testing if SDxSC9 Interaction Improves Model Fit

I calculated this based on this MPlus page: https://www.statmodel.com/chidiff.shtml

https://stats.idre.ucla.edu/mplus/faq/how-can-i-compute-a-chi-square-test-for-nested-models-with-the-mlr-or-mlm-estimators/#:~:text=Typically%20a%20chi%2Dsquare%20difference,freedom%20between%20the%20two%20models.

It references this paper: Satorra, A., & Bentler, P. M. (2010). Ensuring Positiveness of the Scaled Difference Chi-square Test Statistic. Psychometrika, 75(2), 243–248. https://doi.org/10.1007/s11336-009-9135-y

To Compute the Satorra-Bentler Scaled Chi-Square Difference Test:

To calculate the scaling factor:

cd=(p0 * c0 - p1*c1)/(p0 - p1) Where p=parameters and c=correction factor

Then, the test statistic is as follows TRd = -2*(L0 - L1)/cd Where L=loglikelihood

X0 values are the ones with fewer freely estimated parameters X1 values are for the moderation model

This is compared to the chi-square distribution based on the degrees of freedom difference between the two models, in this case, it is 2. Therefore, our model with the interactions fits better than the main effects models.

Compared to ME Only Model

p0=97
p1=99

c0=1.0271
c1=1.0268

cd = (p0*c0-p1*c1)/(p0-p1)
cd
## [1] 1.01225
L0=-50172.965
L1=-50169.544

TRd = -2*(L0-L1)/cd
TRd
## [1] 6.7592
df_paf=p1-p0
qchisq(.95, df_paf)
## [1] 5.991465
pchisq(TRd, df=df_paf, lower.tail=FALSE)
## [1] 0.03406108

Compared to Model Where SDxSC9 Int. Set to 0

p0=98
p1=99

c0=1.0270
c1=1.0268


cd = (p0*c0-p1*c1)/(p0-p1)
cd
## [1] 1.0072
L0=-50172.381
L1=-50169.544


TRd = -2*(L0-L1)/cd
TRd
## [1] 5.633439
df_paf=p1-p0
qchisq(.95, df_paf)
## [1] 3.841459
pchisq(TRd, df=df_paf, lower.tail=FALSE)
## [1] 0.01762103

Who does this interaction not apply to?

According to the Johnson-Neyman intervals, when SC9 <-1.429 stdev and when social deprivation was > 3.067 stdev, they were outside of the range of the interaction. Let’s see who that applies to.

PAF_mod_data %>% 
  filter(SC9_scaled < -1.270) %>% 
  summarize(n(),
            mean(DepComp_scaled, na.rm = TRUE), 
            min(DepComp_scaled), 
            max(DepComp_scaled), 
            mean(PAF_scaled, na.rm = TRUE), 
            min(PAF_scaled),
            max(PAF_scaled))
## # A tibble: 1 x 7
##   `n()` `mean(DepComp_s… `min(DepComp_sc… `max(DepComp_sc… `mean(PAF_scale…
##   <int>            <dbl>            <dbl>            <dbl>            <dbl>
## 1   364           0.0926            -1.71             4.68           -0.342
## # … with 2 more variables: `min(PAF_scaled)` <dbl>, `max(PAF_scaled)` <dbl>
PAF_mod_data %>% 
  filter(DepComp_scaled > 2.766) %>%  
  summarize(n(),
            mean(SC9_scaled, na.rm = TRUE), 
            min(SC9_scaled, na.rm = TRUE), 
            max(SC9_scaled, na.rm=TRUE), 
            mean(PAF_scaled, na.rm = TRUE), 
            min(PAF_scaled, na.rm = TRUE),
            max(PAF_scaled, na.rm = TRUE))
## # A tibble: 1 x 7
##   `n()` `mean(SC9_scale… `min(SC9_scaled… `max(SC9_scaled… `mean(PAF_scale…
##   <int>            <dbl>            <dbl>            <dbl>            <dbl>
## 1    43           -0.140            -3.15             1.29           -0.322
## # … with 2 more variables: `min(PAF_scaled, na.rm = TRUE)` <dbl>,
## #   `max(PAF_scaled, na.rm = TRUE)` <dbl>
PAF_mod_data %>% 
  filter(DepComp_scaled > 2.776 & SC9_scaled < -1.270) %>%  
  summarize(n(),
            mean(PAF_scaled, na.rm = TRUE), 
            min(PAF_scaled, na.rm = TRUE),
            max(PAF_scaled, na.rm = TRUE))
## # A tibble: 1 x 4
##   `n()` `mean(PAF_scaled, na.rm… `min(PAF_scaled, na.rm… `max(PAF_scaled, na.rm…
##   <int>                    <dbl>                   <dbl>                   <dbl>
## 1     3                   -0.954                   -1.51                  0.0494

Checking model assumptions

(Approximate since model is done in Mplus)

# Checking some model assumptions
mean(sc9xpaf_mod$residuals)
## [1] -2.142609e-17
t.test(sc9xpaf_mod$residuals)
## 
##  One Sample t-test
## 
## data:  sc9xpaf_mod$residuals
## t = -1.2625e-15, df = 3245, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.03327605  0.03327605
## sample estimates:
##     mean of x 
## -2.142609e-17
plot(sc9xpaf_mod)

sred = studres(sc9xpaf_mod)

hist(sred, freq=FALSE, 
     main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40) 
yfit<-dnorm(xfit) 
lines(xfit, yfit)

Robustness Check of the Interaction

We tested if the method for extracting the factor scores made a difference in the interpretation of the interaction. The analyses below used factor scores extracted from CFAs containing ONLY one latent variable at a time. I’m using this as the robustness check rather than the preferred method because I think the overall latent variable moderation model accounts for correlations between the latent variables and these are picked up when extracting factor scores from CFAs with all of the latent variables included in the moderation model.

Indiv_FS_All_Variables = Indiv_FS_All_Variables %>% 
  mutate(PAF_scaled = scale(PAF),
         DepComp_scaled = scale(DepComp),
         SC9_scaled = scale(SC9))

sc9xpaf_mod_indiv = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)

summary(sc9xpaf_mod_indiv)
## 
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4193 -0.6713 -0.0051  0.7105  2.1402 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -0.0006067  0.0182564  -0.033   0.9735    
## DepComp_scaled            -0.1107125  0.0187281  -5.912 3.79e-09 ***
## SC9_scaled                 0.1220655  0.0182720   6.680 2.85e-11 ***
## DepComp_scaled:SC9_scaled -0.0436768  0.0194845  -2.242   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9798 on 2887 degrees of freedom
##   (355 observations deleted due to missingness)
## Multiple R-squared:  0.02972,    Adjusted R-squared:  0.02871 
## F-statistic: 29.47 on 3 and 2887 DF,  p-value: < 2.2e-16
johnson_neyman(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is OUTSIDE the interval [-20.12, -1.21], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-2.89, 1.19]

johnson_neyman(sc9xpaf_mod_indiv, pred = SC9_scaled, modx = DepComp_scaled)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is OUTSIDE the interval [1.35, 22.44], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-1.91, 7.54]

sim_slopes(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, data = Indiv_FS_All_Variables, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is OUTSIDE the interval [-20.124, -1.210], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-2.894, 1.193]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of DepComp_scaled when SC9_scaled = -0.999 (- 1 SD): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.067   0.026   -2.536   0.011
## 
## Slope of DepComp_scaled when SC9_scaled =  0.001 (Mean): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.111   0.019   -5.913   0.000
## 
## Slope of DepComp_scaled when SC9_scaled =  1.000 (+ 1 SD): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.154   0.028   -5.602   0.000
interact_plot(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = Indiv_FS_All_Variables) + theme_classic() + theme(text=element_text(size = 16, family="serif"))

Probing SC9xSD –> Ext9 w.stratification

Ext9_mod_data = Ext9_mod_data %>% 
  mutate(DepComp_scaled = scale(DepComp),
         SC9_scaled = scale(SC9),
         Ext9_scaled = scale(Ext9))


sc9xext9_mod = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data)
sc9ext9_ros = data.frame(xmin = -2.9,
                       xmax = 3.221,
                       ymin = -Inf, 
                       ymax = Inf)

interact_plot(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "Social Deprivation", y.label = "Externalizing Symptoms Age 9", legend.main = "School Connectedness", data = Ext9_mod_data) + theme_classic()  + theme(text=element_text(size = 20, family="serif")) + geom_rect(data=sc9ext9_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.3, inherit.aes=FALSE)  + geom_vline(xintercept = 0, colour="black", linetype = "longdash") 

interact_plot(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, interval = TRUE, x.label = "School Connectedness", y.label = "Externalizing Sx Age 9", legend.main = "Social Deprivation", data = Ext9_mod_data) + theme_classic() 

johnson_neyman(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is INSIDE the interval [-3.325, 106.882], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-3.300, 1.595]

johnson_neyman(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is INSIDE the interval [-105.418, 3.221], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]

sim_slopes(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, data = Ext9_mod_data, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is INSIDE the interval [-3.325, 106.882], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-3.300, 1.595]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of DepComp_scaled when SC9_scaled = -1.000 (- 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.194   0.024    8.191   0.000
## 
## Slope of DepComp_scaled when SC9_scaled = -0.000 (Mean): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.227   0.017   13.617   0.000
## 
## Slope of DepComp_scaled when SC9_scaled =  1.000 (+ 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.259   0.025   10.424   0.000
sim_slopes(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, data = Ext9_mod_data, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is INSIDE the interval [-105.418, 3.221], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of SC9_scaled when DepComp_scaled = -1.000 (- 1 SD): 
## 
##     Est.    S.E.    t val.       p
## -------- ------- --------- -------
##   -0.254   0.024   -10.529   0.000
## 
## Slope of SC9_scaled when DepComp_scaled = -0.000 (Mean): 
## 
##     Est.    S.E.    t val.       p
## -------- ------- --------- -------
##   -0.222   0.017   -13.331   0.000
## 
## Slope of SC9_scaled when DepComp_scaled =  1.000 (+ 1 SD): 
## 
##     Est.    S.E.   t val.       p
## -------- ------- -------- -------
##   -0.189   0.024   -7.746   0.000

Getting coefficients, variances, and covariances for and running the Preacher interaction tool.

vcov(sc9xext9_mod)
##                            (Intercept) DepComp_scaled   SC9_scaled
## (Intercept)               2.767544e-04   1.144093e-06 2.565935e-07
## DepComp_scaled            1.144093e-06   2.772599e-04 2.227941e-05
## SC9_scaled                2.565935e-07   2.227941e-05 2.766394e-04
## DepComp_scaled:SC9_scaled 2.504654e-05   1.430275e-05 3.207774e-06
##                           DepComp_scaled:SC9_scaled
## (Intercept)                            2.504654e-05
## DepComp_scaled                         1.430275e-05
## SC9_scaled                             3.207774e-06
## DepComp_scaled:SC9_scaled              3.131165e-04
summary(sc9xext9_mod)
## 
## Call:
## lm(formula = Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5730 -0.6417 -0.0311  0.5608  5.8561 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                0.002607   0.016636   0.157   0.8755    
## DepComp_scaled             0.226741   0.016651  13.617   <2e-16 ***
## SC9_scaled                -0.221720   0.016632 -13.331   <2e-16 ***
## DepComp_scaled:SC9_scaled  0.032589   0.017695   1.842   0.0656 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9444 on 3242 degrees of freedom
## Multiple R-squared:  0.109,  Adjusted R-squared:  0.1082 
## F-statistic: 132.2 on 3 and 3242 DF,  p-value: < 2.2e-16
psych::describe(Ext9_mod_data$SC9_scaled)
##    vars    n mean sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 3246    0  1   0.09    0.08 1.04 -3.3 1.59  4.89 -0.64     0.14 0.02
psych::describe(Ext9_mod_data$DepComp_scaled)
##    vars    n mean sd median trimmed  mad   min  max range skew kurtosis   se
## X1    1 3246    0  1  -0.15   -0.09 0.89 -2.77 7.54 10.31 1.34      4.3 0.02

Preacher Output - SC9 as Moderator

Your Input

X1 = -2.77 X2 = 7.54 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = 0.002607 X Slope = 0.226741 Z Slope = -0.22172 XZ Slope = 0.032589 df = 3240 alpha = 0.05

Asymptotic (Co)variances

var(b0) 0.00027675 var(b1) 0.00027726 var(b2) 0.00027664 var(b3) 0.00031312 cov(b2,b0) 2.6e-7 cov(b3,b1) 0.0000143

Region of Significance

Z at lower bound of region = -3.3254 Z at upper bound of region = 106.8632 (simple slopes are significant inside this region.)

Simple Intercepts and Slopes at Conditional Values of Z

At Z = cv1… simple intercept = 0.2243(0.0235), t=9.5404, p=0 simple slope = 0.1942(0.0237), t=8.1915, p=0 At Z = cv2… simple intercept = 0.0026(0.0166), t=0.1567, p=0.8755 simple slope = 0.2267(0.0167), t=13.6172, p=0 At Z = cv3… simple intercept = -0.2191(0.0235), t=-9.31, p=0 simple slope = 0.2593(0.0249), t=10.4235, p=0

Simple Intercepts and Slopes at Region Boundaries

Lower Bound…
simple intercept = 0.7399(0.0577), t=12.814, p=0 simple slope = 0.1184(0.0604), t=1.9607, p=0.05 Upper Bound…
simple intercept = -23.6911(1.7775), t=-13.3284, p=0 simple slope = 3.7093(1.8918), t=1.9607, p=0.05

Points to Plot

Line for cv1: From {X=-2.77, Y=-0.3135} to {X=7.54, Y=1.6882} Line for cv2: From {X=-2.77, Y=-0.6255} to {X=7.54, Y=1.7122} Line for cv3: From {X=-2.77, Y=-0.9375} to {X=7.54, Y=1.7362}

Preacher Simple Slopes Plot

xx <- c(-2.77,7.54)   #  <-- change to alter plot dims
yy <- c(-1.4722,1.7362)   #  <-- change to alter plot dims
leg <- c(-2.77,-1.0711)   #  <-- change to alter legend location
x <- c(-2.77,7.54)   #  <-- x-coords for lines
y1 <- c(-0.3135,1.6882)
y2 <- c(-0.6255,1.7122)
y3 <- c(-0.9375,1.7362)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='X',ylab='Y',main='MLR 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

Testing if the SDxSC9 Interaction Improves Fit

Compared to ME Only Model

p0=187
p1=191
c0=1.0505
c1=1.0345

cd = (p0*c0-p1*c1)/(p0-p1)
cd
## [1] 0.2865
L0=-70459.229
L1=-70454.076

TRd = -2*(L0-L1)/cd
TRd
## [1] 35.97208
df_paf=p1-p0
qchisq(.95, df_paf)
## [1] 9.487729
pchisq(TRd, df=df_paf, lower.tail=FALSE)
## [1] 2.932224e-07

Compared to Model Where SDxSC9 Int. Set to 0

p0=190
p1=191
c0=1.0386
c1=1.0345

cd = (p0*c0-p1*c1)/(p0-p1)
cd
## [1] 0.2555
L0=-70456.106
L1=-70454.076

TRd = -2*(L0-L1)/cd
TRd
## [1] 15.89041
df_paf=p1-p0
qchisq(.95, df_paf)
## [1] 3.841459
pchisq(TRd, df=df_paf, lower.tail=FALSE)
## [1] 6.711791e-05

Who does this interaction not apply to?

Ext9_mod_data %>% 
   filter(DepComp_scaled > 2.493) %>%  
   summarize(n(),
             mean(SC9_scaled, na.rm = TRUE), 
             min(SC9_scaled, na.rm = TRUE), 
             max(SC9_scaled, na.rm=TRUE), 
             mean(Ext9_scaled, na.rm = TRUE), 
             min(Ext9_scaled, na.rm = TRUE),
             max(Ext9_scaled, na.rm = TRUE))
## # A tibble: 1 x 7
##   `n()` `mean(SC9_scale… `min(SC9_scaled… `max(SC9_scaled… `mean(Ext9_scal…
##   <int>            <dbl>            <dbl>            <dbl>            <dbl>
## 1    58           -0.191            -3.14             1.28            0.506
## # … with 2 more variables: `min(Ext9_scaled, na.rm = TRUE)` <dbl>,
## #   `max(Ext9_scaled, na.rm = TRUE)` <dbl>
Ext9_mod_data %>% 
   filter(SC9_scaled < -2.749) %>%  
   summarize(n(),
             mean(DepComp_scaled, na.rm = TRUE), 
             min(DepComp_scaled, na.rm = TRUE), 
             max(DepComp_scaled, na.rm=TRUE), 
             mean(Ext9_scaled, na.rm = TRUE), 
             min(Ext9_scaled, na.rm = TRUE),
             max(Ext9_scaled, na.rm = TRUE))
## # A tibble: 1 x 7
##   `n()` `mean(DepComp_s… `min(DepComp_sc… `max(DepComp_sc… `mean(Ext9_scal…
##   <int>            <dbl>            <dbl>            <dbl>            <dbl>
## 1    40           0.0132            -1.49             2.81            0.657
## # … with 2 more variables: `min(Ext9_scaled, na.rm = TRUE)` <dbl>,
## #   `max(Ext9_scaled, na.rm = TRUE)` <dbl>
Ext9_mod_data %>% 
  filter(DepComp_scaled > 2.493 & SC9_scaled < -2.749) %>%  
  summarize(n(),
            mean(Ext9_scaled, na.rm = TRUE), 
            min(Ext9_scaled, na.rm = TRUE),
            max(Ext9_scaled, na.rm = TRUE))
## # A tibble: 1 x 4
##   `n()` `mean(Ext9_scaled, na.r… `min(Ext9_scaled, na.r… `max(Ext9_scaled, na.r…
##   <int>                    <dbl>                   <dbl>                   <dbl>
## 1     1                    0.801                   0.801                   0.801

Checking model assumptions

(Approximate since model is done in Mplus)

# Checking some model assumptions
mean(sc9xext9_mod$residuals)
## [1] -5.269587e-18
t.test(sc9xext9_mod$residuals)
## 
##  One Sample t-test
## 
## data:  sc9xext9_mod$residuals
## t = -3.1806e-16, df = 3245, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.0324847  0.0324847
## sample estimates:
##     mean of x 
## -5.269587e-18
plot(sc9xext9_mod)

sred = studres(sc9xext9_mod)

hist(sred, freq=FALSE, 
     main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40) 
yfit<-dnorm(xfit) 
lines(xfit, yfit)

Robustness Check – Alternative Way of Extracting Factor Scores

Indiv_FS_All_Variables = Indiv_FS_All_Variables %>% 
  mutate(Ext9_scaled = scale(EXT9))

sc9xext9_mod_indiv = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)

johnson_neyman(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is INSIDE the interval [-3.27, 99.13], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-2.89, 1.19]

johnson_neyman(sc9xext9_mod_indiv, pred = SC9_scaled, modx = DepComp_scaled)
## JOHNSON-NEYMAN INTERVAL 
## 
## When DepComp_scaled is INSIDE the interval [-56.59, 1.77], the slope of
## SC9_scaled is p < .05.
## 
## Note: The range of observed values of DepComp_scaled is [-1.91, 7.54]

sim_slopes(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, data = Indiv_FS_All_Variables, digits = 3)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SC9_scaled is INSIDE the interval [-3.273, 99.125], the slope of
## DepComp_scaled is p < .05.
## 
## Note: The range of observed values of SC9_scaled is [-2.894, 1.193]
## 
## SIMPLE SLOPES ANALYSIS 
## 
## Slope of DepComp_scaled when SC9_scaled = -0.999 (- 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.207   0.026    7.928   0.000
## 
## Slope of DepComp_scaled when SC9_scaled =  0.002 (Mean): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.242   0.018   13.114   0.000
## 
## Slope of DepComp_scaled when SC9_scaled =  1.002 (+ 1 SD): 
## 
##    Est.    S.E.   t val.       p
## ------- ------- -------- -------
##   0.278   0.027   10.219   0.000
interact_plot(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Externalizing Symptoms (Age 9)", legend.main = "School Connectedness", data = Indiv_FS_All_Variables) + theme_classic() + theme(text=element_text(size = 16, family="serif"))

Hand plotting interaction

PAF_mod_data = PAF_mod_data %>% 
  mutate(SC_bin = case_when(
    SC9_scaled <= -1 ~ 'Low',
    SC9_scaled <= 1 ~ 'Med',
    SC9_scaled > 1 ~ 'High',
    TRUE ~ 'NA'
  ))

Ext9_mod_data = Ext9_mod_data %>% 
  mutate(SC_bin = case_when(
    SC9_scaled <= -1 ~ 'Low',
    SC9_scaled <= 1 ~ 'Med',
    SC9_scaled > 1 ~ 'High',
    TRUE ~ 'NA'
  ))


PAF_mod_data_bin = PAF_mod_data %>% 
  filter(SC_bin != 'NA')

Ext9_mod_data_bin = Ext9_mod_data %>% 
    filter(SC_bin != 'NA')
PAF_mod_data_bin %>%
  ggplot(aes(x=DepComp_scaled, y=PAF_scaled, group=SC_bin, color=SC_bin)) +
  stat_smooth(se = TRUE) +
  xlab("Social Deprivation (Centered)") + ylab("Positive Function (Centered)")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

PAF_mod_data_bin %>%
  ggplot(aes(x=DepComp_scaled, y=PAF_scaled, group=SC_bin, color=SC_bin)) +
  stat_smooth(method = "lm", se = TRUE) +
  xlab("Social Deprivation (Centered)") + ylab("Positive Function (Centered)")
## `geom_smooth()` using formula 'y ~ x'

interact_mod_PAF = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data_bin)

interact_plot(interact_mod_PAF, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = PAF_mod_data_bin) + theme_classic() + theme(text=element_text(size = 16, family="serif"))

Ext9_mod_data_bin %>%
  ggplot(aes(x=DepComp_scaled, y=Ext9_scaled, group=SC_bin, color=SC_bin)) +
  stat_smooth(se = TRUE) +
  xlab("Social Deprivation (Centered)") + ylab("Externalizing Age 9 (Centered)")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Ext9_mod_data_bin %>%
  ggplot(aes(x=DepComp_scaled, y=Ext9_scaled, group=SC_bin, color=SC_bin)) +
  stat_smooth(method="lm",se = TRUE) +
  xlab("Social Deprivation (Centered)") + ylab("Externalizing Age 9 (Centered)")
## `geom_smooth()` using formula 'y ~ x'

interact_mod_Ext9 = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data_bin)

interact_plot(interact_mod_Ext9, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation (Centered)", y.label = "Externalizing Age 9 (Centered)", legend.main = "School Connectedness", data = Ext9_mod_data_bin) + theme_classic() + theme(text=element_text(size = 16, family="serif"))

Who is the Primary Caregiver?

Age 3: IH3int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather

Age5: IH5int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather

Age 9: pcg5idstat 61 Bio Mother
63 Other
62 Bio Father

PC_Info = read_csv('primary_caregiver_relationships.csv')
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   ff_id = col_double(),
##   IH3int5 = col_double(),
##   IH5int5 = col_double(),
##   pcg5idstat = col_double()
## )
PC_Info %>% 
  group_by(IH3int5) %>% 
  summarize(n())
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 7 x 2
##   IH3int5 `n()`
##     <dbl> <int>
## 1       1  3191
## 2       2    15
## 3       3    10
## 4       4     3
## 5       5     4
## 6       6     5
## 7      NA  1670
PC_Info %>% 
  group_by(IH5int5) %>% 
  summarize(n())
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 7 x 2
##   IH5int5 `n()`
##     <dbl> <int>
## 1       1  2914
## 2       2    27
## 3       3    19
## 4       4     7
## 5       5    10
## 6       6     1
## 7      NA  1920
PC_Info %>% 
  group_by(pcg5idstat) %>% 
  summarize(n())
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 4 x 2
##   pcg5idstat `n()`
##        <dbl> <int>
## 1         61  3828
## 2         62   180
## 3         63   194
## 4         NA   696

Hypothetical Resilience Graphs

Adapted from Luthar et al., 2000

Risk = c(2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4)
Competence = c(2,1,3,2,4,3,
               2,1,3,1.5,5,2.5,
               4,2,4.5,3,5,5, 
               2,1,3,3,4,5)
Moderator = c("Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High")
Type = c("promotive","promotive","promotive","promotive","promotive","promotive","prot-reac","prot-reac","prot-reac","prot-reac","prot-reac","prot-reac","prot-stab","prot-stab","prot-stab","prot-stab","prot-stab","prot-stab", "prot-enh", "prot-enh", "prot-enh", "prot-enh", "prot-enh", "prot-enh")

hypothetical_plots = data.frame(Risk, Competence, Moderator, Type)

hypothetical_plots$Moderator <- factor(hypothetical_plots$Moderator, levels = c("Low","Medium","High"))
P = hypothetical_plots %>% 
  filter(Type == "promotive") %>% 
  ggplot(aes(Risk, Competence, group=Moderator)) +
  geom_line(aes(linetype=Moderator), color="black") +
  scale_x_continuous(limits = c(1.9,4.1)) +
  scale_y_continuous(limits = c(0.9,4.1)) +
  ggtitle("Promotive") +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face="bold"), legend.position="none", text=element_text(family="serif"))

PE = hypothetical_plots %>% 
  filter(Type == "prot-enh") %>% 
  ggplot(aes(Risk, Competence, group=Moderator)) +
  geom_line(aes(linetype=Moderator), color="black") +
  scale_x_continuous(limits = c(1.9,4.1)) +
  scale_y_continuous(limits = c(0.9,5.1)) +
  ggtitle("Protective-Enhancing") +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), legend.position="none", text=element_text(family="serif"))

PR = hypothetical_plots %>% 
  filter(Type == "prot-reac") %>% 
  ggplot(aes(Risk, Competence, group=Moderator)) +
  geom_line(aes(linetype=Moderator), color="black") +
  scale_x_continuous(limits = c(1.9,4.1)) +
  scale_y_continuous(limits = c(0.9,5.1)) +
  ggtitle("Protective-Reactive") +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), legend.position="none", text=element_text(family="serif"))

PS = hypothetical_plots %>% 
  filter(Type == "prot-stab") %>% 
  ggplot(aes(Risk, Competence, group=Moderator)) +
  geom_line(aes(linetype=Moderator), color="black") +
  scale_x_continuous(limits = c(1.9,4.1)) +
  scale_y_continuous(limits = c(0.9,5.1)) +
  ggtitle("Protective-Stabilizing") +
  theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), text=element_text(family="serif"), legend.text=element_text(size = 12), legend.position = "bottom", legend.title = element_text(size=12, face="bold"), legend.key = element_rect(fill = "white")) +
  guides(linetype = guide_legend(title.position="top", title.hjust = 0.5), override.aes = list())
Hyp_figure <- ggarrange(P, PE, PR, PS, common.legend = TRUE,legend.grob = get_legend(PS), legend="bottom",
                    ncol = 2, nrow = 2)

Hyp_figure